LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Cleanup code in progress.

Browse Source
This commit is contained in:
salman 2013-02-20 15:02:02 +00:00
parent abca9a346e
commit 4362477a1f
1 changed files with 146 additions and 165 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

311
trimsp/src/trimsp7l.F
View File

@ -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