diff --git a/trimsp/src/trimsp7l.F b/trimsp/src/trimsp7l.F index 1ef405a..2cce871 100644 --- a/trimsp/src/trimsp7l.F +++ b/trimsp/src/trimsp7l.F @@ -1331,7 +1331,6 @@ C C INCREMENT OF DAMAGE, CASCADE AND PHONON ENERGY C DO 70 IV=1,IH1 -C IF(X(IV).LT.0.OR.X(IV).GT.TT) GO TO 70 I=MAX0(MIN0(IDINT(X1(IV)/CW+1.D0),100),1) DENT(I)=DENT(I)+DENS(IV) DMGN(I)=DMGN(I)+DEN(IV) @@ -1349,13 +1348,13 @@ C IF(X(IV).LT.0.OR.X(IV).GT.TT) GO TO 70 GO TO 70 28 PHON(I)=PHON(I)+DEN(IV) ELP(I,JJJ(IV))=ELP(I,JJJ(IV))+DEN(IV) - 70 CONTINUE - DO 80 IV=1,IH1 + 70 CONTINUE + DO IV=1,IH1 ICDI=ICDI+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.DI(JJJ(IV)))) ICSUMS=ICSUMS+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.SB(1))) ICSUM=ICSUM+IDINT(CVMGT(1.D0,0.D0,DENS(IV).GT.0.D0)) - 80 CONTINUE - DO 72 IV=1,IH1 + ENDDO + DO IV=1,IH1 DEN2=DEN(IV)*DEN(IV) DEN3=DEN2*DEN(IV) EEL=EEL+DEN(IV) @@ -1380,9 +1379,9 @@ C IF(X(IV).LT.0.OR.X(IV).GT.TT) GO TO 70 EPL6=EPL6+CVMGT(DEN3*DEN3,0.D0,DEN(IV).LT.DI(JJJ(IV))) ENUCL(IV)=ENUCL(IV)+DENS(IV) EINEL(IV)=EINEL(IV)+DEE(IV) - 72 CONTINUE + ENDDO IF(KK0.EQ.0) GO TO 89 - DO 71 IV=1,IH1 + DO IV=1,IH1 DEWC=DENS(IV)-DEN(IV) DEWC2=DEWC*DEWC DEWC3=DEWC2*DEWC @@ -1392,10 +1391,9 @@ C IF(X(IV).LT.0.OR.X(IV).GT.TT) GO TO 70 EELWC4=EELWC4+DEWC2*DEWC2 EELWC5=EELWC5+DEWC3*DEWC2 EELWC6=EELWC6+DEWC3*DEWC3 - 71 CONTINUE - 89 CONTINUE -C -C IF IRL=0 NO RECOILS ARE FOLLOWED + ENDDO + 89 CONTINUE + IF(IRL.EQ.0) GO TO 27 C C VECTORIZED RECOIL LOOP @@ -1405,38 +1403,31 @@ C C PRIMARY KNOCK-ON ATOMS C DO 6 IV=1,IH1 -cc IF(DEN(IV).LE.SFE) GO TO 6 - IF(DEN(IV).LE.ERC) GO TO 6 - IF(X1(IV).GT.RD.AND.X1(IV).LT.RT) GO TO 6 -C -C CALL NEWREC(NREC1,DEN(IV),X(IV),Y(IV),Z(IV), -C 1 CX(IV),CY(IV),CZ(IV),SX(IV), -C 2 CT(IV),ST(IV),PHI(IV),P(IV), -C 3 ER(1,1),XR(1,1),YR(1,1),ZR(1,1),PHIR(1,1),PSIR(1,1), -C 4 CSXR(1,1),CSYR(1,1),CSZR(1,1),SNXR(1,1),L(1,1) - NREC1=NREC1+1 - ER(NREC1,1)=DEN(IV)-EP(JJJ(IV)) - XR(NREC1,1)=X1(IV) - YR(NREC1,1)=Y(IV)-P(IV)*(SPHI(IV)*CZ(IV) - * -CPHI(IV)*CY(IV)*CX(IV))/SX(IV) - ZR(NREC1,1)=Z(IV)+P(IV)*(SPHI(IV)*CY(IV) - * +CPHI(IV)*CX(IV)*CZ(IV))/SX(IV) - CSXR(NREC1,1)=CX(IV) - CSYR(NREC1,1)=CY(IV) - CSZR(NREC1,1)=CZ(IV) - SNXR(NREC1,1)=SX(IV) - CPHIR(NREC1,1)=-CPHI(IV) - SPHIR(NREC1,1)=-SPHI(IV) - CT(IV)=DMIN1(CT(IV),.99999999D0) - TAR=ST(IV)/(1.D0-CT(IV)) - TAR2=1./DSQRT(1.D0+TAR*TAR) - CPSIR(NREC1,1)=TAR2 - SPSIR(NREC1,1)=TAR*TAR2 - TAUPSR(NREC1,1)=0.D0 - JJR(NREC1,1)=JJJ(IV) - KO(NREC1,JJR(NREC1,1),1)=1 - INOUT(NREC1,1)=SIGN(1.D0,CX(IV)) - NPA=NPA+1 + IF(DEN(IV).LE.ERC) GO TO 6 + IF(X1(IV).GT.RD.AND.X1(IV).LT.RT) GO TO 6 + NREC1=NREC1+1 + ER(NREC1,1)=DEN(IV)-EP(JJJ(IV)) + XR(NREC1,1)=X1(IV) + YR(NREC1,1)=Y(IV)-P(IV)*(SPHI(IV)*CZ(IV)-CPHI(IV)*CY(IV)*CX(IV) + & )/SX(IV) + ZR(NREC1,1)=Z(IV)+P(IV)*(SPHI(IV)*CY(IV)+CPHI(IV)*CX(IV)*CZ(IV) + & )/SX(IV) + CSXR(NREC1,1)=CX(IV) + CSYR(NREC1,1)=CY(IV) + CSZR(NREC1,1)=CZ(IV) + SNXR(NREC1,1)=SX(IV) + CPHIR(NREC1,1)=-CPHI(IV) + SPHIR(NREC1,1)=-SPHI(IV) + CT(IV)=DMIN1(CT(IV),.99999999D0) + TAR=ST(IV)/(1.D0-CT(IV)) + TAR2=1./DSQRT(1.D0+TAR*TAR) + CPSIR(NREC1,1)=TAR2 + SPSIR(NREC1,1)=TAR*TAR2 + TAUPSR(NREC1,1)=0.D0 + JJR(NREC1,1)=JJJ(IV) + KO(NREC1,JJR(NREC1,1),1)=1 + INOUT(NREC1,1)=SIGN(1.D0,CX(IV)) + NPA=NPA+1 6 CONTINUE C IF(NREC1.LT.NUM) GO TO 27 @@ -1449,23 +1440,22 @@ C 1 ,CPSIR(1,1),SPSIR(1,1),CPHIR(1,1),SPHIR(1,1),NREC1) C C MOVE TARGET RECOIL ATOMS TO LIST 2 -CDIR$ IVDEP - DO 91 IREC1=1,NREC1 - IREC=IREC1+NREC2 - ER(IREC,2)=ER(IREC1,1) - XR(IREC,2)=XR(IREC1,1) - YR(IREC,2)=YR(IREC1,1) - ZR(IREC,2)=ZR(IREC1,1) - CSXR(IREC,2)=CSXR(IREC1,1) - CSYR(IREC,2)=CSYR(IREC1,1) - CSZR(IREC,2)=CSZR(IREC1,1) - SNXR(IREC,2)=SNXR(IREC1,1) - TAUPSR(IREC,2)=TAUPSR(IREC1,1) - CPSIR(IREC,2)=CPSIR(IREC1,1) - JJR(IREC,2)=JJR(IREC1,1) - KO(IREC,JJR(IREC,2),2)=KO(IREC1,JJR(IREC1,1),1) - INOUT(IREC,2)=INOUT(IREC1,1) - 91 CONTINUE + DO IREC1=1,NREC1 + IREC=IREC1+NREC2 + ER(IREC,2)=ER(IREC1,1) + XR(IREC,2)=XR(IREC1,1) + YR(IREC,2)=YR(IREC1,1) + ZR(IREC,2)=ZR(IREC1,1) + CSXR(IREC,2)=CSXR(IREC1,1) + CSYR(IREC,2)=CSYR(IREC1,1) + CSZR(IREC,2)=CSZR(IREC1,1) + SNXR(IREC,2)=SNXR(IREC1,1) + TAUPSR(IREC,2)=TAUPSR(IREC1,1) + CPSIR(IREC,2)=CPSIR(IREC1,1) + JJR(IREC,2)=JJR(IREC1,1) + KO(IREC,JJR(IREC,2),2)=KO(IREC1,JJR(IREC1,1),1) + INOUT(IREC,2)=INOUT(IREC1,1) + ENDDO C NREC2=NREC2+NREC1 MAXA=MAX0(MAXA,NREC2) @@ -1474,74 +1464,62 @@ C IF(NREC2.GT.2000) GO TO 8885 GO TO 8886 8885 WRITE(6,8887) - 8887 FORMAT(1X,'ERROR : DIMENSION IN THE RECOIL LOOP ', - 1 'MUST BE INCREASED') -cTR 8887 FORMAT(1X,'ERROR : DIMENSION IN THE RECOIL LOOP MUST BE -cTR 1 INCREASED') + 8887 FORMAT(1X,'ERROR : DIMENSION IN THE RECOIL LOOP ', + & 'MUST BE INCREASED') 8886 CONTINUE C C PROCESS THE PARTICLES IN LIST 2 C C FIND LAYER C - DO 68 IREC1=1,NREC2 - LRR(IREC1,2)=MIN0(ISRCHFGT(L,XX(1),1,XR(IREC1,2)),L) - 68 CONTINUE + DO IREC1=1,NREC2 + LRR(IREC1,2)=MIN0(ISRCHFGT(L,XX(1),1,XR(IREC1,2)),L) + ENDDO C C MOVE PARTICLES C - DO 60 IREC1=1,NREC2 - XR(IREC1,2)=XR(IREC1,2)+(LM(LRR(IREC1,2)) - * +TAUPSR(IREC1,2))*CSXR(IREC1,2) - YR(IREC1,2)=YR(IREC1,2)+(LM(LRR(IREC1,2)) - * +TAUPSR(IREC1,2))*CSYR(IREC1,2) - ZR(IREC1,2)=ZR(IREC1,2)+(LM(LRR(IREC1,2)) - * +TAUPSR(IREC1,2))*CSZR(IREC1,2) - 60 CONTINUE -C - DO 81 IREC1=1,NREC2 - CXR(IREC1)=CSXR(IREC1,2) - CYR(IREC1)=CSYR(IREC1,2) - CZR(IREC1)=CSZR(IREC1,2) - SXR(IREC1)=SNXR(IREC1,2) - DEERS(IREC1)=0.D0 - TS(IREC1)=0.D0 - 81 CONTINUE -C + DO IREC1=1,NREC2 + XR(IREC1,2)=XR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2)) + & *CSXR(IREC1,2) + YR(IREC1,2)=YR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2)) + & *CSYR(IREC1,2) + ZR(IREC1,2)=ZR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2)) + & *CSZR(IREC1,2) + ENDDO + + DO IREC1=1,NREC2 + CXR(IREC1)=CSXR(IREC1,2) + CYR(IREC1)=CSYR(IREC1,2) + CZR(IREC1)=CSZR(IREC1,2) + SXR(IREC1)=SNXR(IREC1,2) + DEERS(IREC1)=0.D0 + TS(IREC1)=0.D0 + ENDDO + KK2=KK0R DO 235 KKR=KK2,0,-1 C C CHOICE OF COLLISION PARTNERS C - DO 303 IREC1=1,NREC2 - call ranlux(random, 1) - JJR(IREC1,1) = ISRCHFGE(NJ(LRR(IREC1,2)),COM(1,LRR(IREC1,2)) -CC 1 ,1,RANF())+JT(LRR(IREC1,2)) -CC 1 ,1,DRAND48())+JT(LRR(IREC1,2)) -CC 1 ,1,DBLE(RAN(ISEED)))+JT(LRR(IREC1,2)) - 1 ,1,DBLE(random))+JT(LRR(IREC1,2)) - 303 CONTINUE -C -CDIR$ IVDEP - DO 236 IREC1=1,NREC2 -CC PHIPR=PI2*RANF() -CC PHIPR=PI2*DRAND48() -CC PHIPR=PI2*DBLE(RAN(ISEED)) - call ranlux(ran2, 2) - PHIPR=PI2*DBLE(ran2(1)) - CPHIR(IREC1,2)=DCOS(PHIPR) - SPHIR(IREC1,2)=DSIN(PHIPR) -CC PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(RANF()+KKR) -CC PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DRAND48()+KKR) -CC PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DBLE(RAN(ISEED))+KKR) - PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DBLE(ran2(2))+KKR) - X2(IREC1)=XR(IREC1,2)-PR(IREC1)*CPHIR(IREC1,2)*SXR(IREC1) - PR(IREC1)=CVMGT(1.D10,PR(IREC1),X2(IREC1).LT.0.0D0 - 1 .OR.X2(IREC1).GT.TT) - BR(IREC1)=PR(IREC1)/A(JJR(IREC1,2),JJR(IREC1,1)) - EPSR(IREC1)=F(JJR(IREC1,2),JJR(IREC1,1))*ER(IREC1,2) - 236 CONTINUE -C + DO IREC1=1,NREC2 + call ranlux(random, 1) + JJR(IREC1,1) = ISRCHFGE(NJ(LRR(IREC1,2)),COM(1,LRR(IREC1,2)) + & ,1,DBLE(random))+JT(LRR(IREC1,2)) + ENDDO + + DO IREC1=1,NREC2 + call ranlux(ran2, 2) + PHIPR=PI2*DBLE(ran2(1)) + CPHIR(IREC1,2)=DCOS(PHIPR) + SPHIR(IREC1,2)=DSIN(PHIPR) + PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DBLE(ran2(2))+KKR) + X2(IREC1)=XR(IREC1,2)-PR(IREC1)*CPHIR(IREC1,2)*SXR(IREC1) + PR(IREC1)=CVMGT(1.D10,PR(IREC1),X2(IREC1).LT.0.0D0 .OR.X2(IREC1 + & ).GT.TT) + BR(IREC1)=PR(IREC1)/A(JJR(IREC1,2),JJR(IREC1,1)) + EPSR(IREC1)=F(JJR(IREC1,2),JJR(IREC1,1))*ER(IREC1,2) + ENDDO + CALL SCOPY(NREC2,BR,1,RR,1) IVMIN=1 IVMAX=NREC2 @@ -1551,67 +1529,70 @@ C IF(IPOTR.NE.1) GO TO 4201 C KR-C POTENTIAL C CALL MAGICKRC(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) - 205 DO 206 IV=IVMIN,IVMAX - IF(RR(IV).LT.1.D-20)THEN - WRITE(99,'(A70)')'in DO 205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - RR(IV)=0.00001D0 - ENDIF - EX1R(IV)=DEXP(-.278544D0*RR(IV)) - EX2R(IV)=DEXP(-.637174D0*RR(IV)) - EX3R(IV)=DEXP(-1.919249D0*RR(IV)) - RRR1=1./RR(IV) - VR(IV)=(.190945D0*EX1R(IV)+.473674D0*EX2R(IV) - # +.335381D0*EX3R(IV))*RRR1 - FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) - V1R(IV)=-(VR(IV)+.053186584080D0*EX1R(IV) - # +.301812757276D0*EX2R(IV)+.643679648869D0*EX3R(IV))*RRR1 - FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ - 1 EPSR(IV)-1.D0 - Q=FR/FR1 - RR(IV)=RR(IV)-Q - TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 - 206 CONTINUE + 205 DO IV=IVMIN,IVMAX + IF(RR(IV).LT.1.D-20)THEN + WRITE(99,'(A70)' + & )'in DO 205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + RR(IV)=0.00001D0 + ENDIF + EX1R(IV)=DEXP(-.278544D0*RR(IV)) + EX2R(IV)=DEXP(-.637174D0*RR(IV)) + EX3R(IV)=DEXP(-1.919249D0*RR(IV)) + RRR1=1./RR(IV) + VR(IV)=(.190945D0*EX1R(IV)+.473674D0*EX2R(IV)+.335381D0 + & *EX3R(IV))*RRR1 + FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) + V1R(IV)=-(VR(IV)+.053186584080D0*EX1R(IV)+.301812757276D0 + & *EX2R(IV)+.643679648869D0*EX3R(IV))*RRR1 + FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/EPSR(IV) + & -1.D0 + Q=FR/FR1 + RR(IV)=RR(IV)-Q + TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 + ENDDO C GET MAX AND MIN INDEX OF TEST FAILURES IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) IF(IVMIN.GT.IVMAX) GO TO 207 IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) IF(IVMIN.GT.IVMAX) GO TO 207 GO TO 205 - 207 DO 208 IV=1,NREC2 - ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) - SQE=DSQRT(EPSR(IV)) - CC=(.235800D0+SQE)/(.126000D0+SQE) - AA=2.D0*EPSR(IV)*(1.D0+(1.0144D0/SQE))*BR(IV)**CC - FF=(DSQRT(AA*AA+1.D0)-AA)*((69350.D0+EPSR(IV)) - # /(83550.D0+EPSR(IV))) - DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) - C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) - C2R(IV)=DMIN1(1.0D0,C*C) - 208 S2R(IV)=1.D0-C2R(IV) + 207 DO IV=1,NREC2 + ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) + SQE=DSQRT(EPSR(IV)) + CC=(.235800D0+SQE)/(.126000D0+SQE) + AA=2.D0*EPSR(IV)*(1.D0+(1.0144D0/SQE))*BR(IV)**CC + FF=(DSQRT(AA*AA+1.D0)-AA)*((69350.D0+EPSR(IV)) /(83550.D0 + & +EPSR(IV))) + DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) + C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) + C2R(IV)=DMIN1(1.0D0,C*C) + S2R(IV)=1.D0-C2R(IV) + ENDDO GO TO 4203 -C + 4201 IF(IPOTR.NE.2) GO TO 4202 C MOLIERE POTENTIAL C CALL MAGICMOL(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) - 4205 DO 4206 IV=IVMIN,IVMAX - IF(RR(IV).LT.1.D-20)THEN - WRITE(99,'(A70)')'in DO 4205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - RR(IV)=0.00001D0 - ENDIF - EX1R(IV)=DEXP(-.3D0*RR(IV)) - EX2R(IV)=DEXP(-1.2D0*RR(IV)) - EX3R(IV)=DEXP(-6.0D0*RR(IV)) - RRR1=1.D0/RR(IV) - VR(IV)=(.35D0*EX1R(IV)+.55D0*EX2R(IV)+.10D0*EX3R(IV))*RRR1 - FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) - V1R(IV)=-(VR(IV)+.105D0*EX1R(IV)+ - # .66D0*EX2R(IV)+.6D0*EX3R(IV))*RRR1 - FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ - 1 EPSR(IV)-1.D0 - Q=FR/FR1 - RR(IV)=RR(IV)-Q - TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 - 4206 CONTINUE + 4205 DO IV=IVMIN,IVMAX + IF(RR(IV).LT.1.D-20)THEN + WRITE(99,'(A70)' + & )'in DO 4205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + RR(IV)=0.00001D0 + ENDIF + EX1R(IV)=DEXP(-.3D0*RR(IV)) + EX2R(IV)=DEXP(-1.2D0*RR(IV)) + EX3R(IV)=DEXP(-6.0D0*RR(IV)) + RRR1=1.D0/RR(IV) + VR(IV)=(.35D0*EX1R(IV)+.55D0*EX2R(IV)+.10D0*EX3R(IV))*RRR1 + FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) + V1R(IV)=-(VR(IV)+.105D0*EX1R(IV)+ .66D0*EX2R(IV)+.6D0*EX3R(IV)) + & *RRR1 + FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ EPSR(IV) + & -1.D0 + Q=FR/FR1 + RR(IV)=RR(IV)-Q + TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 + ENDDO C GET MAX AND MIN INDEX OF TEST FAILURES IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) IF(IVMIN.GT.IVMAX) GO TO 4207