PROGRAM ARCULV C$STDUNIT C The line above was introduced by josef fritscher to port this C program to the Sequent Balance C C AR ION ON CU TARGET WITH SURFACE 100 C MAXIMUM NUMBER OF PARTICLES 2000, C LAST PARTICLE IS ION WHIH ORIGINATES C IN THE X-Y PLANE WITH Z=0 C C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XN(10000),VN(10000),XN1(10000),VN1(10000) DIMENSION AB(10000),AN(10000),AN1(10000),VTABA1(603) DIMENSION FTABA1(603),VTABA2(303),XL(10000),EXTIM(3000) DIMENSION FTABA2(303),VTABA3(303),FTABA3(303) DIMENSION VTABB1(1203),FTABB1(1203),MOVING(3000) DIMENSION XARRAY(9),YARRAY(9),VORZ(9) C**************************************************************** C VTAB(602) TABLE FOR POTENTIAL C FTAB(602) TABLE FOR FORCE C LX(6000) LATTICE LOCATION STORAGE C XN(6000) PARTICLE LOCATION AT STEP N C XN1(6000) PARTICLE LOCATION AT STEP N+1 C VN(6000) PARTICLE VELOCITIES AT STEP N C VN1(6000) PARTICLE VELOCITIES AT STEP N+1 C AB(6000) PARTICLE ACCELERATION AT STEP N-1 C AN(6000) PARTICLE ACCELERATION AT STEP N C AN1(6000) PARTICLE ACCELERATION AT STEP N+1 C N NUMBER OF PARTICLES INCLUDING ION C M N*3 C MM MM-3 C PMA PARTICLE MASS OF TARGETATOMS C PMI PARTICLE MASS OF ION C DE DISSOCIATION ENERGY (EV) C BETA PARAMETER OF MORSE POTENTIAL C RE EQUILIBRIUM DISTANCE IN A C FMIN MINIMUM ENERGY TO ACTIVATE ATOM, GET IT INTO THE C MOVING LIST C DISTANCE 2*C*DISTANCE =MAXIMUM DISTANCE A PARTICLE MOVES C IN A TIMESTEP C H TIMESTEP DELTA T C RMAXSQ MAXIMUM POTENTIAL RADIUS **2 (= CC**2) !!!!!!!! C CC CUT OF RADIUS =C*FCUTOF C FCUTOF CUT-OFF FACTOR C C 1/2 LATTICE CONSTANT = 1.8075 FOR CU C FMOV MOVING FACTOR C CMOV MOVING RADIUS C MCENT M-INDEX AND C NXCENT*C COORDINATES OF SURFACE CENTER PARTICLE C NYCENT*C FOR WHICH THE SURFACE BINDING ENERGY C NZCENT*C IS CALCULATED C FX ACCELERATION PARTICLE EXPERIENCES IN X DIRECTION C FY ACCELERATION PARTICLE EXPERIENCES IN Y DIRECTION C FZ ACCELERATION PARTICLE EXPERIENCES IN Z DIRECTION C TMAX MAXIMUM TIME OF CALCULATION C PRINTN NUMBER OF STORED TIMESTEPS C PM DE BETA RE C FMIN DISTANCE TMAX PRINTN C UNITS....UNITS.....UNITS C XN/XN1 ANGSTROM=A C VN/VN1 ANGSTROM/FSEC C AB/AN/AN1 ANGSTROM/FSEC**2 C ENERGY EV C FX,FY,FZ ANGSTROM/FSEC**2 C****************************************************************** C C FORMATTABLE C C****************************************************************** 1000 FORMAT(I8) 1058 FORMAT(G15.7) 1004 FORMAT(1X,4G15.7) 1077 FORMAT(3G15.7) 1078 FORMAT(1X,I8,3G15.7) 1159 FORMAT(1X,4E15.7) 1052 FORMAT(4G15.7) 1987 FORMAT(I8,4G15.7) 1008 FORMAT(1X,3F15.7,/,1X,3F15.7) C******************************************************************* C C READ INPUT DATA, CONSTANTS, FACTORS USED ARE DEFINED C C******************************************************************* WRITE(6,1001) 1001 FORMAT(1X,' NX=; NUMBER OF PARTICLE IN X DIRECTION') READ(5,*) NXOR WRITE(6,1056) 1056 FORMAT(1X,' NY=; NUMBER OF PARTICLE IN Y DIRECTION') READ(5,*) NYOR WRITE(6,1057) 1057 FORMAT(1X,' NZ=; NUMBER OF PARTICLE IN Z DIRECTION') READ(5,*) NZOR WRITE(6,9702) 9702 FORMAT(1X,' READ: EION= ION ENERGY IN EV',/) READ(5,*) EION WRITE(6,1160) 1160 FORMAT(1X,/,' READ: FMIN(EV),DISTANCE,TMAX,PRINTN') READ(5,*) FMIN,DISTANCE,TMAX,PRINTN WRITE(6,4567) 4567 FORMAT(1X,' TAKE AWAY CENTER ATOM =1') WRITE(6,4569) 4569 FORMAT(1X,' TAKE AWAY 0.5 MONOLAYER = 10') WRITE(6,4570) 4570 FORMAT(1X,' ATTENTION: CHOOSE NYOR : EVEN') READ(5,*)ITAWAY WRITE(6,201) 201 FORMAT(1X,' NUMBER OF ARRAYS: 1,2,3,-4-,5,6,7,8,-9-') READ(5,*) IARRAY WRITE(6,202) 202 FORMAT(1X,' NUMBER OF POINTS = NUMBER**2') READ(5,*)NUMBER C*********************************************************** C TARGETCONSTANTS: C 0.5 * LATTICE CONSTANT C********************************************************** C=1.8075 XARRAY(1)=0. XARRAY(2)=C/2. XARRAY(3)=C/2. XARRAY(4)=C XARRAY(5)=C XARRAY(6)=C XARRAY(7)=1.5*C XARRAY(8)=1.5*C XARRAY(9)=2.*C YARRAY(1)=0. YARRAY(2)=C/2 YARRAY(3)=C/2 YARRAY(4)=0. YARRAY(5)=C YARRAY(6)=C YARRAY(7)=C/2 YARRAY(8)=C/2 YARRAY(9)=0. VORZ(1)=1. VORZ(2)=1. VORZ(3)=-1. VORZ(4)=1. VORZ(5)=-1. VORZ(6)=1. VORZ(7)=1. VORZ(8)=-1. VORZ(9)=1. C************************************************************** C C TARGETCONSTANTS: C 0.5 * LATTICE CONSTANT C PMI ION MASS-AR C PMA ATOM MASS-CU C EION ION ENERGY C FMIN MINIMUM ACCELERATION AN ATOM C HAS TO EXPERIENCE TO BE SET INTO C MOTION, I.E. PUT ON THE MOVING C PARTICLE LIST C =0.0002 TYPICAL VALUE C 2.*C*DISTANCE MAXIMUM DISTANCE THE FASTEST ATOM C IS ALLOWED TO MOVE DURING A TIMESTEP C =0.015-0.02 TYPICAL VALUE C C*************************************************************** C C=1.8075 CND=1.414213*C PMI=40. PMI2=52.2*PMI PMA=64. PMA2=52.2*PMA PMR=PMA/PMI C C************************************************************** C C OPEN OUTPUTFILE OUTPUT.DAT WRITE DATA TO IT C **************************************************************** C OPEN(UNIT=1,FILE='OUTPUTp.DAT',STATUS='NEW') WRITE(1,9701) 9701 FORMAT(1X,' AR--CU, CU SURFACE (100) 1 MAX.NUMBER OF ATOMS=2000',/) WRITE(1,8853) 8853 FORMAT(1X,/,' ANGLE OF INCIDENCE (WITH TARGET 1 NORMAL) === 0 GRAD ',/) WRITE(1,1913) WRITE(1,*)NXOR 1913 FORMAT(1X,' NUMBER OF PARTICLES IN X DIRECTION= ') WRITE(1,1914) 1914 FORMAT(1X,' NUMBER OF PARTICLES IN Y DIRECTION= ') WRITE(1,*)NYOR WRITE(1,1915) 1915 FORMAT(1X,' NUMBER OF PARTICLES IN Z DIRECTION= ') WRITE(1,*)NZOR WRITE(1,9716) 9716 FORMAT(1X,'AR ION ON CU TARGET WITH SURFACE 100',/, 1 ' EION = ') WRITE(1,*)EION WRITE(1,1006) 1006 FORMAT(1X,/,' FMIN ',8X,'DISTANCE',5X,' TMAX',5X, 1 ' PRINTN') WRITE(1,*)FMIN,DISTANCE,TMAX,PRINTN WRITE(1,9706) 9706 FORMAT(1X,' 0.5*LATTICE UNIT= ', 1 ' CLOSEST NEIGHBOR DISTANCE=',/) WRITE(1,*)C,CND WRITE(1,9201) 9201 FORMAT(1X,' NUMBER OF ARRAYS= ') WRITE(1,*)IARRAY WRITE(1,9202) 9202 FORMAT(1X,' NUMBER OF POINTS/ARRAY=NUMBER**2, NUMBER=') WRITE(1,*)NUMBER C***************************************************************** C POTENTIAL (VTAB) AND FORCE (FTAB) TABLE ARE SET UP C POTENTIAL.DAT ===3 IS THE FILE TO STORE POTENTIAL DATA C C**************************************************************** C POTENTIAL A ATOM-ATOM ATOM-ATOM ATOM - ATOM C C POTENTIAL A1 BORN MAYER A1 BORN MAYER C C V=AA1*EXP(-BA1*RIJ) C AA1=22564EV BA1=5.088/A RCA1=1.5A C POTENTIAL/FORCE TABLE: VTABAI,FTABAI C---------------------------------------------------------------- C OPEN(UNIT=3,FILE='POTENTIAL.DAT',STATUS='NEW') AA1=22565. BA1=-5.088 RCA1=1.5 F2A1=(-1.)*(AA1*BA1*9.58E-3)/PMA RCA1SQ=RCA1**2 KA1MAX=600 DSA1=RCA1SQ/REAL(KA1MAX) KA1SAF=602 DO 2401 I=1,KA1SAF S=REAL(I)*DSA1 RIJ=SQRT(S) VA1=EXP(BA1*RIJ) VTABA1(I)=AA1*VA1 FTABA1(I)=(F2A1/RIJ)*VA1 C WRITE(3,1077)RIJ,VTABA1(I),FTABA1(I) 2401 CONTINUE WRITE(1,9708)AA1,BA1,RCA1 9708 FORMAT(1X,/,' BM-POT. A-A: AA1=',G15.7,'BA1=',G15.7, 1 'CUTOFFRAD=',G15.7) C---------------------------------------------------------------- C POTENTIAL A ATOM-ATOM ATOM-ATOM ATOM - ATOM C C POTENTIAL A2 SPLINE A2 SPLINE C V=CA0+CA1*RIJ+CA2*RIJ**2+CA3*RIJ**3 C CA0=570.5EV CA1=-853.0EV/A CA2=428.0EV/A** CA3=-72.0EV/A**3 C RCA2=1.988A C POTENTIAL/FORCE TABLE: VTABA2,FTABA2 C----------------------------------------------------------------- CA0=570.5 CA1=-853.0 CA2=428.0 CA3=-72.0 RCA2=1.988 F2A2=-9.58E-3/PMA RCA2SQ=RCA2**2 KA2MAX=300 DSA2=RCA2SQ/REAL(KA2MAX) KA2SAF=302 DO 2402 I=1,KA2SAF S=REAL(I)*DSA2 RIJ=SQRT(S) VTABA2(I)=CA0+CA1*RIJ+CA2*RIJ**2+CA3*RIJ**3 FTABA2(I)=(F2A2/RIJ)*(CA1+2.*CA2*RIJ+3.*CA3*RIJ**2) C WRITE(3,1077)RIJ,VTABA2(I),FTABA2(I) 2402 CONTINUE WRITE(1,9709)CA0,CA1,CA2,CA3,RCA2 9709 FORMAT(1X,/,' SPLINE A-A: CA0=',G15.7,'CA1=',G15.7, 1 /,'CA2=',G15.7 1 ' CA3=',G15.7,'CUTOFFRAD=',G15.7) C------------------------------------------------------------------ C POTENTIAL A ATOM-ATOM ATOM-ATOM ATOM - ATOM C C POTENTIAL A3 MORSE POT. MORSE POT. C V=DEA*(EXP(-BETAA2*(RIJ-REA))-2.*EXP(-BETAA*(RIJ-REA))) C DEA=0.48EV BETAA=1.405/A BETAA2=2.*BETAA REA=2.628A C RCA3=4.338A=2.4*C C POTENTIAL/FORCE TABLE: VTABA3,FTABA3 C------------------------------------------------------------------ DEA=0.48 BETAA=1.405 BETAA2=2.*BETAA REA=2.628 RCA3=4.33 F2A3=(BETAA2*DEA*9.58E-3)/PMA RCA3SQ=RCA3**2 KA3MAX=300 DSA3=RCA3SQ/REAL(KA3MAX) KA3SAF=302 DO 2403 I=1,KA3SAF S=REAL(I)*DSA3 RIJ=SQRT(S) VTABA3(I)=DEA*(EXP(-BETAA2*(RIJ-REA))-2.*EXP(-BETAA*(RIJ-REA))) FTABA3(I)=(F2A3/RIJ)*EXP(-BETAA2*(RIJ-REA)) FTABA3(I)=FTABA3(I)-(F2A3/RIJ)*EXP(-BETAA*(RIJ-REA)) C WRITE(3,1077)RIJ,VTABA3(I),FTABA3(I) 2403 CONTINUE VERRA=DEA*(EXP(-BETAA2*(RCA3-REA))-2.*EXP(-BETAA*(RCA3-REA))) WRITE(1,9710) DEA,BETAA,REA,RCA3,VERRA 9710 FORMAT(1X,/,' MORSE P. A-A: DE=',G15.7,'BETA=',G15.7,/, 1 'RE=,'G15.7, 6 'CUTOFFRAD=',G15.7, ' POT AT CUTOFFRAD=',G15.7) C------------------------------------------------------------------- C POTENTIAL B ION-ATOM ION-ATOM ION - ATOM C C POTENTIAL B1 BORN MAYER B1 BORN MAYER C C V=AB*EXP(-BB*RIJ) C AA1=59874EV BB1=7.2/A RCB1=3.0 C POTENTIAL/FORCE TABLE: VTABB1,FTABB1 C C ATTENTION: FTABB1 IS THE ACCELERATION A TARGET ATOM EXPEREIENCES C DUE TO THE ION, TO OBTAIN THE ACCELARATION OF THE ION MULTIPLY C BY PMR: C FTABB1(ION)=FTABB1(I)*PMA/PMI = FTABB1(I)*PMR C C------------------------------------------------------------------- AB1=59874. BB1=-7.2 RCB1=3.00 F2B1=(-1.)*(AB1*BB1*9.58E-3)/PMA RCB1SQ=RCB1**2 KB1MAX=1200 DSB1=RCB1SQ/REAL(KB1MAX) KB1SAF=1202 DO 2404 I=1,KB1SAF S=REAL(I)*DSB1 RIJ=SQRT(S) VB1=EXP(BB1*RIJ) VTABB1(I)=AB1*VB1 FTABB1(I)=(F2B1/RIJ)*VB1 C WRITE(3,1077)RIJ,VTABB1(I),FTABB1(I) 2404 CONTINUE VERRI=AB1*EXP(BB1*RCB1) WRITE(1,9711)AB1,BB1,RCB1,VERRI 9711 FORMAT(1X,/,' BM-POT. I-A: AB1=',G15.7,'BB1=',G15.7, 1 'CUTOFFRAD=',G15.7,/, 6 'POT. AT CUTOFFRAD=',G15.7,/) C CLOSE(UNIT=3) C C********************************************************************* C C LATTICE CONSTRUCTION C C********************************************************************* C NXCENT=NXOR+1 NYCENT=NYOR+1 NZCENT=2 NX=2*NXOR NY=2*NYOR NZ=2*NZOR ML=0 NL=0 DO 512 I=2,NZ IF(ITAWAY.EQ.10.AND.I.EQ.2) THEN JLLI=3 JINCR=2 ELSE JLLI=2 JINCR=1 ENDIF DO 513 J=JLLI,NY,JINCR DO 514 K=2,NX L=I+J+K A1=REAL(L/2) A2=REAL(L)/2. A3=DIM(A2,A1) IF(A3.GT.0.1) GOTO 514 IF(I.EQ.NZCENT.AND.J.EQ.NYCENT.AND.K.EQ.NXCENT) THEN MCENT=ML+1 NCENT=NL+1 ENDIF XL(ML+1)=REAL(K)*C XL(ML+2)=REAL(J)*C XL(ML+3)=REAL(I)*C C WRITE(6,1078) N,XL(M+1),XL(M+2),XL(M+3) ML=ML+3 NL=NL+1 C WRITE(6,1060)NL,ML 514 CONTINUE 513 CONTINUE 512 CONTINUE MML=ML-3 WRITE(6,1060)NL,ML WRITE(1,1060) WRITE(1,*)NL,ML 1060 FORMAT(1X,' NL= ',I8,' ML= ',I8) C C------------------------------------------------------------------ C MOVING LOGIC MOVING LOGIC MOVING LOGIC MOVING LOGIC C C KMOVING(I)= 0 NOT MOVING C 1 MOVING, BUT INSIDE CRYSTAL C C ATOM 2 XXXMAX C 4 YYYMAX C 6 ZZZMAX C C ION 10 XXXMAX C 12 YYYMAX C 14 ZZZMAX C------------------------------------------------------------------ C XXMIN=1.5*C XXMAX=REAL(NX+0.5)*C YYMIN=1.5*C YYMAX=REAL(NY+0.5)*C ZZMIN=1.5*C ZZMAX=REAL(NZ+0.5)*C C C**************************************************************** C CALCULATE SURFACE BINDUNG ENERGY C**************************************************************** SBE=0. C WRITE(6,*)XL(MCENT),XL(MCENT+1),XL(MCENT+2) DO 768 J=1,ML,3 IF(J.EQ.MCENT) GOTO 768 XIJ=XL(J)-XL(MCENT) IF(ABS(XIJ).GT.RCA3) GOTO 768 YIJ=XL(J+1)-XL(MCENT+1) IF(ABS(YIJ).GT.RCA3) GOTO 768 ZIJ=XL(J+2)-XL(MCENT+2) IF(ABS(ZIJ).GT.RCA3) GOTO 768 RIJSQ=XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ IF(RIJSQ.GE.RCA3SQ) GOTO 768 IF(RIJSQ.LE.RCA1SQ) THEN SDS=RIJSQ/DSA1 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA1(K) VK1=VTABA1(K+1) VK2=VTABA1(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 SBE=SBE+VIJ ELSE IF(RIJSQ.LE.RCA2) THEN SDS=RIJSQ/DSA2 K=INT(SDS) XI=SDS-REAL(K) VK=VTABA2(K) VK1=VTABA2(K+1) VK2=VTABA2(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 SBE=SBE+VIJ ELSE SDS=RIJSQ/DSA3 K=INT(SDS) XI=SDS-REAL(K) VK=VTABA3(K) VK1=VTABA3(K+1) VK2=VTABA3(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 SBE=SBE+VIJ ENDIF ENDIF 768 CONTINUE WRITE(1,1096) 1096 FORMAT(1X,/,' SURFACE BINDING ENERGY = .... EV') WRITE(1,*)SBE C C**************************************************************** C C NOW TAKE AWAY ATOM IF ITAWAY=1 C**************************************************************** C WRITE(1,1917) 1917 FORMAT(1X,' TAKE AWAY CENTER ATOM =') WRITE(1,*)ITAWAY WRITE(1,3456) 3456 FORMAT(1X, ' MCENT, NCENT') WRITE(1,*)MCENT,NCENT IF(ITAWAY.EQ.1) THEN DO 7111 IK=MCENT,MML 7111 XL(IK)=XL(IK+3) ML=ML-3 NL=NL-1 MML=MML-3 ENDIF C*************************************************************** C ENERGY OF ION AND TOTAL ENERGY OF SYSTEM ARE CALCULATED C C*************************************************************** IION=ML+1 C IBGIN=JSECNDS(0) ENERGY=0. SUBLIM=0. C C DO LOOP HAS OTGHER UPPER LIMITS THAN DO LOOP C FOR ENERGY AT END OF CASCADE AS AT THE START C OF THE CASCADE THE ION IS NOT INTERACTING\ C WITH THE LATTICE SO POTENTIAL (ION)=0 C IN ADDITION NO MOVING ATOMS EXCEPT FOR ION !!! C IION=ML+1 DO 68 J=1,MML,3 JJ=J+3 XHELP1=XL(J) XHELP2=XL(J+1) XHELP3=XL(J+2) X1LL=XHELP1-RCA3 X1UL=XHELP1+RCA3 X2LL=XHELP2-RCA3 X2UL=XHELP2+RCA3 X3LL=XHELP3-RCA3 X3UL=XHELP3+RCA3 DO 69 I=JJ,ML,3 IF(I.EQ.J) GOTO 69 XXN1=XL(I) IF(XXN1.LT.X1LL) GOTO 69 IF(XXN1.GT.X1UL) GOTO 69 XXN2=XL(I+1) IF(XXN2.LT.X2LL) GOTO 69 IF(XXN2.GT.X2UL) GOTO 69 XXN3=XL(I+2) IF(XXN3.LT.X3LL) GOTO 69 IF(XXN3.GT.X3UL) GOTO 69 XIJ=XHELP1-XXN1 YIJ=XHELP2-XXN2 ZIJ=XHELP3-XXN3 RIJSQ=XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ IF(RIJSQ.GE.RCA3SQ) GOTO 69 IF(RIJSQ.LE.RCA1SQ) THEN SDS=RIJSQ/DSA1 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA1(K) VK1=VTABA1(K+1) VK2=VTABA1(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 ENERGY=ENERGY+VIJ-VERRA SUBLIM=SUBLIM+VIJ ELSE IF(RIJSQ.LE.RCA2) THEN SDS=RIJSQ/DSA2 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA2(K) VK1=VTABA2(K+1) VK2=VTABA2(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 ENERGY=ENERGY+VIJ-VERRA SUBLIM=SUBLIM+VIJ ELSE SDS=RIJSQ/DSA3 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA3(K) VK1=VTABA3(K+1) VK2=VTABA3(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 ENERGY=ENERGY+VIJ-VERRA SUBLIM=SUBLIM+VIJ ENDIF ENDIF 69 CONTINUE 68 CONTINUE C C SUBLIMATIONENERGY IS TAKEN HERE BEFORE ION ENERGY IS C ADDEDD TO ENERGY C SUBLIM=(-1.)*SUBLIM/REAL(NL-1) ENERGY=ENERGY+EION WRITE(1,2405)ENERGY,SUBLIM WRITE(6,2405)ENERGY,SUBLIM 2405 FORMAT(1X,'ENERGY= ',E11.4,' HEAT SUBLIMATION=',F10.3) WRITE(1,*)ENERGY,SUBLIM C IFINAL=JSECNDS(IBGIN) C WRITE(6,1820)IFINAL C1820 FORMAT(1X,' TOTAL ENERGY IS CALCULATED, TIME= ',I8) C****************************************************************** C C START THE BIG DO LOOP C C****************************************************************** IHELP1=0 IHELP2=0 IHELP3=-500 LIMIT1=NUMBER*2-1 LIMIT2=NUMBER W2HX=C/(2.*REAL(NUMBER)) C C HERE !!!!!!+++++++++++++++++++++++++++++++++++++++++++ C DO 1886 INDEX0=1,IARRAY DO 1888 INDEX1=1,LIMIT1 LIHELP=2*NUMBER-INDEX1 DO 1877 INDEX2=1,LIMIT2 IF(INDEX1.LE.NUMBER.AND.INDEX2.GT.INDEX1) GOTO 1877 IF(INDEX1.GT.NUMBER.AND.INDEX2.GT.LIHELP) GOTO 1877 C------------------------------------------- C RESET LATTICE C-------------------------------------------- C DO 1918 I= 1,ML 1918 XN(I)=XL(I) M=ML N=NL C C******************************************************************* C C POINT OF IMPACT -ION VELOCITY / ANGLE OF INCIDENCE C******************************************************************* XX=FLOAT(INDEX1)*W2HX I4=INDEX1+INDEX2 I3=MOD(I4,2) IF(I3.EQ.0) THEN YY=FLOAT(INDEX2)*W2HX-W2HX*0.666666 ELSE YY=FLOAT(INDEX2)*W2HX-W2HX*0.333333 ENDIF XN(M+1)=(REAL(NX)/2.+1.)*C+XARRAY(INDEX0)+XX XN(M+2)=(REAL(NY)/2.+1.)*C+YARRAY(INDEX0) 1 +VORZ(INDEX0)*YY XN(M+3)=0. VN(M+1)=0. VN(M+2)=0. IF(EION.GT.0.) THEN VN(M+3)=SQRT(EION/PMI2) ELSE VN(M+3)=-1.*SQRT(-1.*EION/PMI2) ENDIF M=M+3 N=N+1 MM=M-3 WRITE(6,9705)INDEX0,INDEX1,INDEX2 9705 FORMAT(1X,'ARRAY=',I4'INDEX1= ',I8,' INDEX2= ',I8,/) WRITE(6,1075) XN(MM+1),XN(MM+2),XN(MM+3) 1075 FORMAT(1X,' ION: X=',F10.3,' Y=',F10.3,' Z=',F10.3) C WRITE(6,9713) C9713 FORMAT(1X,/,1X,' ION DATA: VX, VY, VZ') C WRITE(6,*)(VN(MM+I),I=1,3) TIME=0. FMINSQ=FMIN*FMIN VMAXSQ=VN(M)**2+VN(M-1)**2+VN(M-2)**2 VMAX=SQRT(VMAXSQ) C IBEGIN=JSECNDS(0) C C OPEN HER FOR PRINTOUTS !!!!!!! ISI=0 C IF(ISI.EQ.1) THEN WRITE(6,1956) 1956 FORMAT(1X,' NO POSITION PRINTOUT TO FILE',/, 1 ' ENERGY CALCULATION ONLY AT START AND END , YES=1',/, 1 ' ENERGY AND POSITIONS ARE STORED PRINTN TIMES, YES=2 ONLY') READ(5,1000)IBB ELSE IBB=1 ENDIF C**************************************************************** C SETTINGS AB(I) = 0. C AN(I) = 0. C AN1(I) = 0. C VN1(I) = VN(I) C XN1(I) = XN(I) C C STARTING DATA STORED ON FILE C C************************************************************** DO 456 I=1,MM 456 VN(I)=0. DO 19 I=1,M AN(I)=0. AN1(I)=0. AB(I)=0. IF(I.GE.MM) GOTO 19 VN1(I)=VN(I) 19 XN1(I)=XN(I) HB=DISTANCE*2.*C/VMAX TPRIN=TMAX/PRINTN TPRINT=TPRIN IF(IBB.EQ.2) THEN ICOUNT=0 DO 1201 I=1,M,3 ICOUNT=ICOUNT+1 C WRITE(1,1987) ICOUNT,TIME,VN(I),VN(I+1),VN(I+2) 1201 WRITE(1,1987) ICOUNT,TIME,XN(I),XN(I+1),XN(I+2) ENDIF C C*************************************************************** C C STEP I C CALCULATE AN(I) C C THIS STEP IS TURNED OFF UNDER THE FOLLOWING ASSUMPTIONS: C 1) TARGET ATOMS ARE INITIALLY AT EQUILIBRIUM POSITIONS C 2)INCOMING PARTICLE IS ORIGINALLY FAR ENOUGH AWAY TO C EXORT NO FORCE ON ANY ATOM AND VICE VERSA C C*************************************************************** C STEP II C CALCULATE XN1(I) C C**************************************************************** IF(FMIN.GT.0.) THEN DO 831 I=1,N EXTIM(I)=0. 831 MOVING(I)=0 MOVING(N)=1 EXTIM(N)=0. ELSE DO 9715 I=1,N 9715 MOVING(I)=1 ENDIF IMAX=N C**************************************************************** C C NEXT DELTA T STEP STARTS C C**************************************************************** 70 CONTINUE H=DISTANCE*2.*C/VMAX RSTEP=H/HB H1=H**2/6. H2=RSTEP*H1 H6=H/6. H4=(3.+RSTEP)*H1 C C IBGIN=JSECNDS(0) C C WRITE(6,*)N,M DO 400 K=1,N IF(MOVING(K).EQ.0) GOTO 400 I=3*(K-1)+1 XN1(I)=XN(I)+VN(I)*H+AN(I)*H4-AB(I)*H2 XN1(I+1)=XN(I+1)+VN(I+1)*H+H4*AN(I+1)-AB(I+1)*H2 XN1(I+2)=XN(I+2)+VN(I+2)*H+H4*AN(I+2)-AB(I+2)*H2 C WRITE(6,*)K,XN1(I+2),XN(I+2) 400 CONTINUE C WRITE(6,*)XN1(M),VN(M),H,XN(M),AN(M),AB(M) C IFINAL=JSECNDS(IBGIN) C WRITE(6,182)IFINAL C182 FORMAT(1X,' STEP II IS FINISHED, TIME= ',I8) C************************************************************* C C SET UP OF NEIGHBORHOOD LIST C C************************************************************** C IBGIN=JSECNDS(0) C DO 2901 ICELL=1,NCELL C JHEAD(ICELL)=0 C2901 CONTINUE C DO 2931 I=1,N C JJ=3*(I-1)+1 C ICELL=1+INT(XN1(JJ))+INT(XN1(JJ+1))*LCUBE+ C 1 INT(XN1(JJ+2))*LCUBSQ C LIST(I)=JHEAD(ICELL) C JHEAD(ICELL)=I C2931 CONTINUE C IFINAL=JSECNDS(IBGIN) C WRITE(6,2904) IFINAL C2904 FORMAT(1X,' END OF NEIGHBOR LIST, TIME= ',I8) C DO 2932 I=1,NCELL C2932 WRITE(6,*) JHEAD(I) C C*************************************************************** C C STEP III C CALCULATE AN1(I) WITH THE NEW SPACE VECTORS XN1(I), C UPDATE MOVING LIST C C*************************************************************** C C IBGIN=JSECNDS(0) C NN1=N-1 DO 500 K=1,NN1 IF(MOVING(K).EQ.0) GOTO 500 I=3*K-2 FX=0. FY=0. FZ=0. XHELP1=XN1(I) XHELP2=XN1(I+1) XHELP3=XN1(I+2) X1LL=XHELP1-RCA3 X1UL=XHELP1+RCA3 X2LL=XHELP2-RCA3 X2UL=XHELP2+RCA3 X3LL=XHELP3-RCA3 X3UL=XHELP3+RCA3 DO 718 JJ=1,N J=3*JJ-2 IF(I.EQ.J) GOTO 718 XN11=XN1(J) IF(XN11.LT.X1LL) GOTO 718 IF(XN11.GT.X1UL) GOTO 718 XN12=XN1(J+1) IF(XN12.LT.X2LL) GOTO 718 IF(XN12.GT.X2UL) GOTO 718 XN13=XN1(J+2) IF(XN13.LT.X3LL) GOTO 718 IF(XN13.GT.X3UL) GOTO 718 XIJ=XN11-XHELP1 YIJ=XN12-XHELP2 ZIJ=XN13-XHELP3 RIJSQ=XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ IF(JJ.EQ.N) THEN IF(RIJSQ.GE.RCB1SQ) GOTO 718 SDS=RIJSQ/DSB1 KK=INT(SDS) IF(KK.EQ.0) GOTO 7777 XI=SDS-REAL(KK) VK=FTABB1(KK) VK1=FTABB1(KK+1) VK2=FTABB1(KK+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) FF=T1+(T2-T1)*XI*0.5 FFX=XIJ*FF FFY=YIJ*FF FFZ=ZIJ*FF FX=FX-FFX FY=FY-FFY FZ=FZ-FFZ C C NO QUESTION AFTER FMIN, AS IT IS THE ION AND C THEREFORE MOVING C ELSE IF(RIJSQ.GE.RCA3SQ) GOTO 718 IF(RIJSQ.LE.RCA1SQ) THEN SDS=RIJSQ/DSA1 KK=INT(SDS) IF(KK.EQ.0) GOTO 7777 XI=SDS-REAL(KK) VK=FTABA1(KK) VK1=FTABA1(KK+1) VK2=FTABA1(KK+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) FF=T1+(T2-T1)*XI*0.5 FFX=XIJ*FF FFY=YIJ*FF FFZ=ZIJ*FF FX=FX-FFX FY=FY-FFY FZ=FZ-FFZ IF(ABS(FFX).LE.FMIN.AND.ABS(FFY).LE.FMIN.AND.ABS(FFZ).LE.FMIN) 1 GOTO 718 IF(MOVING(JJ).EQ.0) THEN MOVING(JJ)=1 EXTIM(JJ)=REAL(IFIX(TIME)) ENDIF C WRITE(6,*) JJ ELSE IF(RIJSQ.LE.RCA2SQ) THEN SDS=RIJSQ/DSA2 KK=INT(SDS) IF(KK.EQ.0) GOTO 7777 XI=SDS-REAL(KK) VK=FTABA2(KK) VK1=FTABA2(KK+1) FF=VK+(VK1-VK)*XI FFX=XIJ*FF FFY=YIJ*FF FFZ=ZIJ*FF FX=FX-FFX FY=FY-FFY FZ=FZ-FFZ IF(ABS(FFX).LE.FMIN.AND.ABS(FFY).LE.FMIN.AND.ABS(FFZ).LE.FMIN) 1 GOTO 718 IF(MOVING(JJ).EQ.0) THEN MOVING(JJ)=1 EXTIM(JJ)=REAL(IFIX(TIME)) ENDIF C WRITE(6,*) JJ ELSE SDS=RIJSQ/DSA3 KK=INT(SDS) IF(KK.EQ.0) GOTO 7777 XI=SDS-REAL(KK) VK=FTABA3(KK) VK1=FTABA3(KK+1) FF=VK+(VK1-VK)*XI FFX=XIJ*FF FFY=YIJ*FF FFZ=ZIJ*FF FX=FX-FFX FY=FY-FFY FZ=FZ-FFZ IF(ABS(FFX).LE.FMIN.AND.ABS(FFY).LE.FMIN.AND.ABS(FFZ).LE.FMIN) 1 GOTO 718 IF(MOVING(JJ).EQ.0)THEN MOVING(JJ)=1 EXTIM(JJ)=REAL(IFIX(TIME)) ENDIF C WRITE(6,*) JJ ENDIF ENDIF ENDIF 718 CONTINUE AN1(I)=FX AN1(I+1)=FY AN1(I+2)=FZ 500 CONTINUE C-------------------- C NOW FORCE ON ION C-------------------- FX=0. FY=0. FZ=0. XHELP1=XN1(M-2) XHELP2=XN1(M-1) XHELP3=XN1(M) X1LL=XHELP1-RCB1 X1UL=XHELP1+RCB1 X2LL=XHELP2-RCB1 X2UL=XHELP2+RCB1 X3LL=XHELP3-RCB1 X3UL=XHELP3+RCB1 DO 9718 JJ=1,NN1 J=3*JJ-2 XN11=XN1(J) IF(XN11.LT.X1LL) GOTO 9718 IF(XN11.GT.X1UL) GOTO 9718 XN12=XN1(J+1) IF(XN12.LT.X2LL) GOTO 9718 IF(XN12.GT.X2UL) GOTO 9718 XN13=XN1(J+2) IF(XN13.LT.X3LL) GOTO 9718 IF(XN13.GT.X3UL) GOTO 9718 XIJ=XN11-XHELP1 YIJ=XN12-XHELP2 ZIJ=XN13-XHELP3 RIJSQ=XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ IF(RIJSQ.GE.RCB1SQ) GOTO 9718 SDS=RIJSQ/DSB1 KK=INT(SDS) IF(KK.EQ.0) GOTO 7777 XI=SDS-REAL(KK) VK=FTABB1(KK) VK1=FTABB1(KK+1) VK2=FTABB1(KK+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) FF=T1+(T2-T1)*XI*0.5 FFX=XIJ*FF FFY=YIJ*FF FFZ=ZIJ*FF FX=FX-FFX*PMR FY=FY-FFY*PMR FZ=FZ-FFZ*PMR IF(ABS(FFX).LE.FMIN.AND.ABS(FFY).LE.FMIN.AND.ABS(FFZ).LE.FMIN) 1 GOTO 9718 IF(MOVING(JJ).EQ.0) THEN MOVING(JJ)=1 EXTIM(JJ)=REAL(IFIX(TIME)) ENDIF C WRITE(6,*) JJ 9718 CONTINUE AN1(M-2)=FX AN1(M-1)=FY AN1(M)=FZ C IFINAL=JSECNDS(IBGIN) C WRITE(6,678)IFINAL,TIME C678 FORMAT(1X,'END STEP III, TIME= ', I8,' FSEC=',G12.4) C C****************************************************************** C C STEP IV C CALCULATE VN1(I) AND RESTORE DATA C C****************************************************************** C C IBGIN=JSECNDS(0) C VMAXSQ=0. H5=H6*(3.+2.*RSTEP)/(1.+RSTEP) H3=(3.+RSTEP)*H6 H7=RSTEP*RSTEP/(1+RSTEP)*H6 DO 800 K=1,N I=3*(K-1)+1 IF(MOVING(K).EQ.0) GO TO 900 VN1I1=VN(I)+H5*AN1(I)+H3*AN(I)-AB(I)*H7 VN1I2=VN(I+1)+H5*AN1(I+1)+H3*AN(I+1)-AB(I+1)*H7 VN1I3=VN(I+2)+H5*AN1(I+2)+H3*AN(I+2)-AB(I+2)*H7 VISQ=VN1I1*VN1I1+VN1I2*VN1I2+VN1I3*VN1I3 IF(VISQ.GT.VMAXSQ) THEN IF(MOVING(K).EQ.1) THEN VMAXSQ=VISQ IMAX=K ENDIF ENDIF C WRITE(6,*)VISQ VN(I)=VN1I1 VN(I+1)=VN1I2 VN(I+2)=VN1I3 AB(I)=AN(I) AN(I)=AN1(I) AB(I+1)=AN(I+1) AB(I+2)=AN(I+2) AN(I+1)=AN1(I+1) AN(I+2)=AN1(I+2) 900 XN(I)=XN1(I) XN(I+1)=XN1(I+1) 800 XN(I+2)=XN1(I+2) HB=H VMAX=SQRT(VMAXSQ) C WRITE(6,*) VN(M),XN(M) C WRITE(6,*) VMAX IF(IMAX.EQ.N) THEN PM2=PMI2 ELSE PM2=PMA2 ENDIF TIME=TIME+H EMAX=PM2*VMAXSQ C IFINAL=JSECNDS(IBEGIN) C WRITE(6,1989) IFINAL,TIME,H,EMAX C1989 FORMAT(1X,' T= ',I8,' FS= ',G12.4,' H= ',G15.7, C 1 ' EMAX= ',G15.7) C*************************************************************** C CHECK FOR DATA STORING AND C INFO ON TERMINAL END OF CALCULATION C C START OF IF THEN LOOP AAAA C C*************************************************************** IF(TIME.GE.TPRINT.OR.EMAX.LT.SUBLIM) THEN IF(IBB.EQ.2) THEN ICOUNT=0 DO 50 I=1,M,3 C WRITE(6,1008)XN(I),XN(I+1),XN(I+2),VN(I),VN(I+1),VN(I+2) ICOUNT=ICOUNT+1 C WRITE(1,1987)ICOUNT,TIME,VN(I),VN(I+1),VN(I+2) 50 WRITE(1,1987)ICOUNT,TIME,XN(I),XN(I+1),XN(I+2) ELSE IF(EMAX.LT.SUBLIM) IBB=3 ENDIF C***************************************************************** C C CALCULATE ENERGY OF SYSTEM AT TIME T C C***************************************************************** ENERGY=0. DO 60 J=1,MM,3 XHELP1=XN(J) XHELP2=XN(J+1) XHELP3=XN(J+2) C--------------------------- IF(IBB.EQ.1) GOTO 561 C--------------------------- X1LL=XHELP1-RCA3 X1UL=XHELP1+RCA3 X2LL=XHELP2-RCA3 X2UL=XHELP2+RCA3 X3LL=XHELP3-RCA3 X3UL=XHELP3+RCA3 JJ=J+3 DO 61 I=JJ,M,3 XXN1=XN(I) IF(XXN1.LT.X1LL) GOTO 61 IF(XXN1.GT.X1UL) GOTO 61 XXN2=XN(I+1) IF(XXN2.LT.X2LL) GOTO 61 IF(XXN2.GT.X2UL) GOTO 61 XXN3=XN(I+2) IF(XXN3.LT.X3LL) GOTO 61 IF(XXN3.GT.X3UL) GOTO 61 XIJ=XHELP1-XXN1 YIJ=XHELP2-XXN2 ZIJ=XHELP3-XXN3 RIJSQ=XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ IF(I.EQ.IION) THEN IF(RIJSQ.GE.RCB1SQ) GOTO 61 SDS=RIJSQ/DSB1 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABB1(K) VK1=VTABB1(K+1) VK2=VTABB1(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 ENERGY=ENERGY+VIJ-VERRI C WRITE(6,*)IION,J,J,I,VIJ,K ELSE IF(RIJSQ.GE.RCA3SQ) GOTO 61 IF(RIJSQ.LE.RCA1SQ) THEN SDS=RIJSQ/DSA1 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA1(K) VK1=VTABA1(K+1) VK2=VTABA1(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 ENERGY=ENERGY+VIJ-VERRA C WRITE(6,*)IION,J,I,I,VIJ,K ELSE IF(RIJSQ.LE.RCA2) THEN SDS=RIJSQ/DSA2 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA2(K) VK1=VTABA2(K+1) VK2=VTABA2(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 ENERGY=ENERGY+VIJ-VERRA C WRITE(6,*)IION,,I,VIJ,VIJ,K ELSE SDS=RIJSQ/DSA3 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA3(K) VK1=VTABA3(K+1) VK2=VTABA3(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 ENERGY=ENERGY+VIJ-VERRA C WRITE(6,*)IION,J,I,VIJ,K,K ENDIF ENDIF ENDIF C WRITE(6,*)RIJSQ,K,VIJ,ENERGY 61 CONTINUE 561 CONTINUE ICHECK=(J+2)/3 IF(MOVING(ICHECK).EQ.1) THEN IF(XHELP1.LT.XXMIN) MOVING(ICHECK)=2 IF(XHELP1.GT.XXMAX) MOVING(ICHECK)=3 IF(XHELP2.LT.YYMIN) MOVING(ICHECK)=4 IF(XHELP2.GT.YYMAX) MOVING(ICHECK)=5 IF(XHELP3.LT.ZZMIN) MOVING(ICHECK)=6 IF(XHELP3.GT.ZZMAX) MOVING(ICHECK)=7 IF(MOVING(ICHECK).GE.2) THEN EXTIM(ICHECK)=EXTIM(ICHECK)+TIME/1000. ENDIF ENDIF CCC IF(MOVING(ICHECK).GE.2) WRITE(6,8014) ICHECK,MOVING(ICHECK), CCC 1 TIME,EMAX 8014 FORMAT(1X,' ICHECK=',2I8,' TIME/ENERGY',2G15.7) ENERGY=ENERGY+PMA2*(VN(J)**2+VN(J+1)**2+VN(J+2)**2) C WRITE(6,*)J,ENERGY,VN(J),VN(J+1),VN(J+2) C WRITE(6,*)VN(J),VN(J+1),VN(J+2),ENERGY 60 CONTINUE ENERGY=ENERGY+PMI2*(VN(MM+1)**2+VN(MM+2)**2+VN(MM+3)**2) IF(MOVING(N).EQ.1) THEN IF(XN(MM+1).LT.XXMIN) MOVING(N)=10 IF(XN(MM+1).GT.XXMAX) MOVING(N)=11 IF(XN(MM+2).LT.YYMIN) MOVING(N)=12 IF(XN(MM+2).GT.YYMAX) MOVING(N)=13 IF(XN(MM+3).LT.ZZMIN.AND.VN(MM+3).LT.0.) MOVING(N)=14 IF(XN(MM+3).GT.ZZMAX) MOVING(N)=15 IF(MOVING(N).GT.1) THEN EXTIM(N)=EXTIM(N)+TIME/1000. ENDIF ENDIF CCC IF(MOVING(N).GE.2) WRITE(6,8014) N,MOVING(N), CCC 1 TIME,EMAX C C**************************************************************** C C IBB =1 SUPRESS PRINTOUT OF ALL ATOMS (X,V) C =2 PRINTOUT OF XN,VN PRINTN TIMES IS DONE C =3 IS SET TO 3 AS SOON AS EMAX < SUBLIM C CAUSES FINAL ENERGY CALCULATION AND PRINTOUT C C*************************************************************** TPRINT=TPRINT+TPRIN ENDIF C----------------------------------------------------------------- C C END OF IF...THEN..LOOP AAAA, RETURN TO 70 JUMP C C------------------------------------------------------------------ C IFINAL=JSECNDS(IBGIN) C WRITE(6,1823)IFINAL C1823 FORMAT(1X,' STEP IV IS FINISHED, TIME= ',I8) C IF(IBB.NE.3)THEN IF(TIME.LT.TMAX) GOTO 70 ELSE WRITE(6,1010) 1010 FORMAT(' RUN FINISHED !!!!! RUN FINISHED!!!!') ENDIF C**************************************************************** C C CALCULATIONN FOR ONE IMPACT IS FINISHED C FINAL PRINTOUT TO 1 C C*************************************************************** IHEL=0 DO 8999 I=1,N IF(MOVING(I).GE.1) IHEL=IHEL+1 IF(MOVING(I).GT.1) THEN BINDING=0. II=3*I-2 XHELP1=XN(II) XHELP2=XN(II+1) XHELP3=XN(II+2) X1LL=XHELP1-RCA3 X1UL=XHELP1+RCA3 X2LL=XHELP2-RCA3 X2UL=XHELP2+RCA3 X3LL=XHELP3-RCA3 X3UL=XHELP3+RCA3 DO 2800 J=1,MM,3 IF(J.EQ.II) GOTO 2800 XXN1=XN(J) IF(XXN1.LT.X1LL) GOTO 2800 IF(XXN1.GT.X1UL) GOTO 2800 XXN2=XN(J+1) IF(XXN2.LT.X2LL) GOTO 2800 IF(XXN2.GT.X2UL) GOTO 2800 XXN3=XN(J+2) IF(XXN3.LT.X3LL) GOTO 2800 IF(XXN3.GT.X3UL) GOTO 2800 XIJ=XXN1-XHELP1 YIJ=XXN2-XHELP2 ZIJ=XXN3-XHELP3 RIJSQ=XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ IF(J.EQ.IION.OR.I.EQ.N) THEN IF(RIJSQ.GE.RCB1SQ) GOTO 2800 SDS=RIJSQ/DSB1 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABB1(K) VK1=VTABB1(K+1) VK2=VTABB1(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 BINDING=BINDING+VIJ C WRITE(6,*)IION,J,J,I,VIJ,K ELSE IF(RIJSQ.GE.RCA3SQ) GOTO 2800 IF(RIJSQ.LE.RCA1SQ) THEN SDS=RIJSQ/DSA1 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA1(K) VK1=VTABA1(K+1) VK2=VTABA1(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 BINDING=BINDING+VIJ ELSE IF(RIJSQ.LE.RCA2) THEN SDS=RIJSQ/DSA2 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA2(K) VK1=VTABA2(K+1) VK2=VTABA2(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 BINDING=BINDING+VIJ ELSE SDS=RIJSQ/DSA3 K=INT(SDS) IF(K.EQ.0) GOTO 7777 XI=SDS-REAL(K) VK=VTABA3(K) VK1=VTABA3(K+1) VK2=VTABA3(K+2) T1=VK+(VK1-VK)*XI T2=VK1+(VK2-VK1)*(XI-1.0) VIJ=T1+(T2-T1)*XI*0.5 BINDING=BINDING+VIJ C WRITE(6,*)IION,J,I,VIJ,K,K ENDIF ENDIF ENDIF C WRITE(6,*)RIJSQ,K,VIJ,BINDING 2800 CONTINUE IF(I.EQ.N) THEN PM2=PMI2 ELSE PM2=PMA2 ENDIF ESPUT=PM2*(VN(II)**2+VN(II+1)**2+VN(II+2)**2) ETOTAL=ESPUT+BINDING IF(ETOTAL.GE.0.) THEN WRITE(1,3567) INDEX0,INDEX1,INDEX2 3567 FORMAT(1X,3I8) WRITE(1,3568) I,MOVING(I),ETOTAL,ESPUT 3568 FORMAT(1X,2I8,2E15.7) WRITE(1,*) (XN(II+K),K=0,2) WRITE(1,*) (VN(II+K),K=0,2) TOTIME=IFIX(EXTIM(I)) EMTIME=(EXTIM(I)-TOTIME)*1000. WRITE(1,*) TOTIME,EMTIME ENDIF ENDIF 8999 CONTINUE WRITE(1,3567)INDEX0,IHELP1,IHELP3 WRITE(1,*) IHEL WRITE(1,3567)INDEX0,IHELP1,IHELP2 WRITE(1,1030) 1030 FORMAT(1X,'ENERGY= TIME= EMAX= ') WRITE(1,*) ENERGY,TIME,EMAX C IFINAL=JSECNDS(IBEGIN) C FINAL=REAL(IFINAL)/60. C WRITE(6,7771)FINAL C7771 FORMAT(1X,' TIME= ',G12.2, 'MINUTES') 7777 WRITE(6,7778) CONTINUE 1877 CONTINUE 1887 CONTINUE 1888 CONTINUE 1886 CONTINUE 7779 WRITE(6,7778) 7778 FORMAT(1X,' HELP HELP IMPACT PARAMETER TO CLOSE FOR FTAB') CLOSE(UNIT=1) STOP END