hdiff output

r30245/djwgr1.f90 2016-04-07 10:30:08.967750859 +0100 r30244/djwgr1.f90 2016-04-07 10:30:09.163753475 +0100
 21: !  NRIGIDBODY = # rigid bodies 21: !  NRIGIDBODY = # rigid bodies
 22: ! 22: !
 23: SUBROUTINE DJWGR1(NATOMS,X,V,ENERGY,GTEST) 23: SUBROUTINE DJWGR1(NATOMS,X,V,ENERGY,GTEST)
 24: USE GENRIGID 24: USE GENRIGID
 25: USE COMMONS, ONLY : NHEXAMERS 25: USE COMMONS, ONLY : NHEXAMERS
 26: IMPLICIT NONE 26: IMPLICIT NONE
 27: LOGICAL GTEST 27: LOGICAL GTEST
 28: INTEGER NATOMS, J1, J2, J3, J4, NPOS1, NPOS2 28: INTEGER NATOMS, J1, J2, J3, J4, NPOS1, NPOS2
 29: DOUBLE PRECISION X(3*NATOMS), V(3*NATOMS), ENERGY, DUMMY2, DUMMY3, DUMMY, DIST, SIGMA, XDUMM, RHO, RDIST, RAD, EPSEFF 29: DOUBLE PRECISION X(3*NATOMS), V(3*NATOMS), ENERGY, DUMMY2, DUMMY3, DUMMY, DIST, SIGMA, XDUMM, RHO, RDIST, RAD, EPSEFF
 30: DOUBLE PRECISION RADHEX, SIGMAHEX, SIGMAPH 30: DOUBLE PRECISION RADHEX, SIGMAHEX, SIGMAPH
 31: DOUBLE PRECISION FATT, DFATT, FREP, DFREP 31: DOUBLE PRECISION FATT, DFATT, FREP, DFREP, FREPHEX, DFREPHEX, FREPPH, DFREPPH
 32: ! 32: !
 33: ! Derivatives of the pairwise site-site terms in terms of distance 33: ! Derivatives of the pairwise site-site terms in terms of distance
 34: ! 34: !
 35: FATT(RHO,XDUMM)=-1.0D0 + (1.0D0 - EXP(RHO*(1.0D0 - XDUMM)))**2 35: FATT(RHO,XDUMM)=-1.0D0 + (1.0D0 - EXP(RHO*(1.0D0 - XDUMM)))**2
 36: DFATT(RHO,XDUMM)=2.0D0*(-EXP(2.0D0*RHO*(1.0D0-XDUMM)) + EXP(RHO*(1.0D0-XDUMM)))*RHO 36: DFATT(RHO,XDUMM)=2.0D0*(-EXP(2.0D0*RHO*(1.0D0-XDUMM)) + EXP(RHO*(1.0D0-XDUMM)))*RHO
 37: FREP(SIGMA,XDUMM)=(SIGMA/XDUMM)**12 37: FREP(SIGMA,XDUMM)=(SIGMA/XDUMM)**12
 38: DFREP(SIGMA,XDUMM)=-12.0D0*(SIGMA/XDUMM)**12/XDUMM 38: DFREP(SIGMA,XDUMM)=-12.0D0*(SIGMA/XDUMM)**12/XDUMM
  39: FREPHEX(SIGMA,XDUMM)=(SIGMAHEX/XDUMM)**12
  40: DFREPHEX(SIGMA,XDUMM)=-12.0D0*(SIGMAHEX/XDUMM)**12/XDUMM
  41: FREPPH(SIGMA,XDUMM)=(SIGMAPH/XDUMM)**12
  42: DFREPPH(SIGMA,XDUMM)=-12.0D0*(SIGMAPH/XDUMM)**12/XDUMM
 39:  43: 
 40: ENERGY=0.0D0 44: ENERGY=0.0D0
 41: IF (GTEST) V(1:3*NATOMS)=0.0D0 45: IF (GTEST) V(1:3*NATOMS)=0.0D0
 42: ! 46: !
 43: ! 5 Morse plus two axial site pentamers from 47: ! 5 Morse plus two axial site pentamers from
 44: ! S.N. Fejer, T. James, J. Hernandez-Rojas and D.J. Wales, Phys. Chem. Chem. Phys., 11, 2098-2104 (2009).  48: ! S.N. Fejer, T. James, J. Hernandez-Rojas and D.J. Wales, Phys. Chem. Chem. Phys., 11, 2098-2104 (2009). 
 45: ! Energy Landscapes for Shells Assembled from Pentagonal and Hexagonal Pyramids  49: ! Energy Landscapes for Shells Assembled from Pentagonal and Hexagonal Pyramids 
 46: ! 50: !
 47: RAD=5.0D0 51: RAD=5.0D0
 48: RADHEX=RAD*2.0*0.5877852522924731D0  ! 2 * Sin[36] to give the same edge length 52: RADHEX=RAD*2.0*0.5877852522924731D0  ! 2 * Sin[36] to give the same edge length
 49:  53: 
 50: RHO=3.0D0 54: RHO=3.0D0
 51: SIGMA=(1.0D0+RAD*SQRT((5.0D0+SQRT(5.0D0))/2.0D0)) 55: SIGMA=(1.0D0+RAD*SQRT((5.0D0+SQRT(5.0D0))/2.0D0))
 52: SIGMAHEX=(1.0D0+RADHEX*SQRT((5.0D0+SQRT(5.0D0))/2.0D0)) 56: SIGMAHEX=(1.0D0+RADHEX*SQRT((5.0D0+SQRT(5.0D0))/2.0D0))
 53: SIGMAPH=0.5D0*(SIGMA + SIGMAHEX) 57: SIGMAPH=0.5D0*(SIGMA + SIGMAHEX)
 54: EPSEFF=0.4D0 58: EPSEFF=0.4D0
  59: 
  60: DO J1=1,NRIGIDBODY
  61:    DO J2=J1+1,NRIGIDBODY
 55: ! 62: !
 56: ! Three different sorts of axial repulsion 63: ! Three different sorts of axial repulsion
  64: ! Same for all hexamers and pentamers.
 57: ! 65: !
 58: ! pent-pent first 
 59: ! 
 60: DO J1=1,NRIGIDBODY-NHEXAMERS 
 61:    DO J2=J1+1,NRIGIDBODY-NHEXAMERS 
 62:       NPOS1=RIGIDGROUPS(1,J1) 66:       NPOS1=RIGIDGROUPS(1,J1)
 63:       NPOS2=RIGIDGROUPS(1,J2) 67:       NPOS2=RIGIDGROUPS(1,J2)
 64:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 68:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2)
 65:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term 69:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term
 66:       IF (GTEST) THEN 70:       IF (GTEST) THEN
  71:          IF (J1.LE.NRIGIDBODY-NHEXAMERS) THEN ! J1 is a pentamer
  72:             IF (J2.LE.NRIGIDBODY-NHEXAMERS) THEN ! J2 is a pentamer
  73:                DUMMY2=EPSEFF*DFREP(SIGMA,DIST)
  74:             ELSE ! J2 is a hexamer
  75:                DUMMY2=EPSEFF*DFREPPH(SIGMAPH,DIST)
  76:             ENDIF
  77:          ELSE ! J2 is a hexamer, and J2 > J1 is also a hexamer
  78:             DUMMY2=EPSEFF*DFREPHEX(SIGMAHEX,DIST)
  79:          ENDIF
 67:          RDIST=1.0D0/DIST 80:          RDIST=1.0D0/DIST
 68:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST)*RDIST 
 69:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 81:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)
 70:       ENDIF 82:       ENDIF
 71:       NPOS1=RIGIDGROUPS(1,J1) 83:       NPOS1=RIGIDGROUPS(1,J1)
 72:       NPOS2=RIGIDGROUPS(2,J2) 84:       NPOS2=RIGIDGROUPS(2,J2)
 73:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 85:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2)
 74:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term 86:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term
 75:       IF (GTEST) THEN 87:       IF (GTEST) THEN
  88:          IF (J1.LE.NRIGIDBODY-NHEXAMERS) THEN ! J1 is a pentamer
  89:             IF (J2.LE.NRIGIDBODY-NHEXAMERS) THEN ! J2 is a pentamer
  90:                DUMMY2=EPSEFF*DFREP(SIGMA,DIST)
  91:             ELSE ! J2 is a hexamer
  92:                DUMMY2=EPSEFF*DFREPPH(SIGMAPH,DIST)
  93:             ENDIF
  94:          ELSE ! J2 is a hexamer, and J2 > J1 is also a hexamer
  95:             DUMMY2=EPSEFF*DFREPHEX(SIGMAHEX,DIST)
  96:          ENDIF
 76:          RDIST=1.0D0/DIST 97:          RDIST=1.0D0/DIST
 77:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST)*RDIST 
 78:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 98:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)
 79:       ENDIF 99:       ENDIF
 80:       NPOS1=RIGIDGROUPS(2,J1)100:       NPOS1=RIGIDGROUPS(2,J1)
 81:       NPOS2=RIGIDGROUPS(1,J2)101:       NPOS2=RIGIDGROUPS(1,J2)
 82:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2)102:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2)
 83:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial 2-axial repulsive term103:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial 2-axial repulsive term
 84:       IF (GTEST) THEN104:       IF (GTEST) THEN
 105:          IF (J1.LE.NRIGIDBODY-NHEXAMERS) THEN ! J1 is a pentamer
 106:             IF (J2.LE.NRIGIDBODY-NHEXAMERS) THEN ! J2 is a pentamer
 107:                DUMMY2=EPSEFF*DFREP(SIGMA,DIST)
 108:             ELSE ! J2 is a hexamer
 109:                DUMMY2=EPSEFF*DFREPPH(SIGMAPH,DIST)
 110:             ENDIF
 111:          ELSE ! J2 is a hexamer, and J2 > J1 is also a hexamer
 112:             DUMMY2=EPSEFF*DFREPHEX(SIGMAHEX,DIST)
 113:          ENDIF
 85:          RDIST=1.0D0/DIST114:          RDIST=1.0D0/DIST
 86:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST)*RDIST 
 87:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 
 88:       ENDIF 
 89:    ENDDO 
 90: ENDDO 
 91: ! 
 92: ! pent-hex second 
 93: ! 
 94: DO J1=1,NRIGIDBODY-NHEXAMERS 
 95:    DO J2=NRIGIDBODY-NHEXAMERS+1,NRIGIDBODY 
 96:       NPOS1=RIGIDGROUPS(1,J1) 
 97:       NPOS2=RIGIDGROUPS(1,J2) 
 98:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 
 99:       ENERGY=ENERGY+EPSEFF*FREP(SIGMAPH,DIST)  ! axial-axial repulsive term 
100:       IF (GTEST) THEN 
101:          RDIST=1.0D0/DIST 
102:          DUMMY2=EPSEFF*DFREP(SIGMAPH,DIST)*RDIST 
103:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 
104:       ENDIF 
105:       NPOS1=RIGIDGROUPS(1,J1) 
106:       NPOS2=RIGIDGROUPS(2,J2) 
107:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 
108:       ENERGY=ENERGY+EPSEFF*FREP(SIGMAPH,DIST)  ! axial-axial repulsive term 
109:       IF (GTEST) THEN 
110:          RDIST=1.0D0/DIST 
111:          DUMMY2=EPSEFF*DFREP(SIGMAPH,DIST)*RDIST 
112:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 
113:       ENDIF 
114:       NPOS1=RIGIDGROUPS(2,J1) 
115:       NPOS2=RIGIDGROUPS(1,J2) 
116:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 
117:       ENERGY=ENERGY+EPSEFF*FREP(SIGMAPH,DIST)  ! axial 2-axial repulsive term 
118:       IF (GTEST) THEN 
119:          RDIST=1.0D0/DIST 
120:          DUMMY2=EPSEFF*DFREP(SIGMAPH,DIST)*RDIST 
121:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 
122:       ENDIF 
123:    ENDDO 
124: ENDDO 
125: ! 
126: ! hex-hex third 
127: ! 
128: DO J1=NRIGIDBODY-NHEXAMERS+1,NRIGIDBODY 
129:    DO J2=J1+1,NRIGIDBODY 
130:       NPOS1=RIGIDGROUPS(1,J1) 
131:       NPOS2=RIGIDGROUPS(1,J2) 
132:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 
133:       ENERGY=ENERGY+EPSEFF*FREP(SIGMAHEX,DIST)  ! axial-axial repulsive term 
134:       IF (GTEST) THEN 
135:          RDIST=1.0D0/DIST 
136:          DUMMY2=EPSEFF*DFREP(SIGMAHEX,DIST)*RDIST 
137:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 
138:       ENDIF 
139:       NPOS1=RIGIDGROUPS(1,J1) 
140:       NPOS2=RIGIDGROUPS(2,J2) 
141:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 
142:       ENERGY=ENERGY+EPSEFF*FREP(SIGMAHEX,DIST)  ! axial-axial repulsive term 
143:       IF (GTEST) THEN 
144:          RDIST=1.0D0/DIST 
145:          DUMMY2=EPSEFF*DFREP(SIGMAHEX,DIST)*RDIST 
146:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V) 
147:       ENDIF 
148:       NPOS1=RIGIDGROUPS(2,J1) 
149:       NPOS2=RIGIDGROUPS(1,J2) 
150:       DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2) 
151:       ENERGY=ENERGY+EPSEFF*FREP(SIGMAHEX,DIST)  ! axial 2-axial repulsive term 
152:       IF (GTEST) THEN 
153:          RDIST=1.0D0/DIST 
154:          DUMMY2=EPSEFF*DFREP(SIGMAHEX,DIST)*RDIST 
155:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)115:          CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)
156:       ENDIF116:       ENDIF
157:    ENDDO 
158: ENDDO 
159: !117: !
160: ! Sum over the attractive sites.118: ! Sum over the attractive sites.
161: ! Same for p-p, p-h, and h-h.119: ! Same for p-p, p-h, and h-h.
162: !120: !
163: DO J1=1,NRIGIDBODY 
164:    DO J2=J1+1,NRIGIDBODY 
165:       DO J3=3,NSITEPERBODY(J1)     ! # Morse sites in rb J1 NSITEPERBODY(J1)121:       DO J3=3,NSITEPERBODY(J1)     ! # Morse sites in rb J1 NSITEPERBODY(J1)
166:          NPOS1=RIGIDGROUPS(J3,J1)    ! where is this site in the list?122:          NPOS1=RIGIDGROUPS(J3,J1)    ! where is this site in the list?
167:          DO J4=3,NSITEPERBODY(J2)  ! # Morse sites in rb J2 NSITEPERBODY(J2)123:          DO J4=3,NSITEPERBODY(J2)  ! # Morse sites in rb J2 NSITEPERBODY(J2)
168:             NPOS2=RIGIDGROUPS(J4,J2) ! where is this site in the list?124:             NPOS2=RIGIDGROUPS(J4,J2) ! where is this site in the list?
169:             DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2)125:             DIST=SQRT((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))**2+(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))**2+(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))**2)
170:             ENERGY=ENERGY+FATT(RHO,DIST) 126:             ENERGY=ENERGY+FATT(RHO,DIST) 
171:             IF (GTEST) THEN127:             IF (GTEST) THEN
 128:                DUMMY2=DFATT(RHO,DIST)
172:                RDIST=1.0D0/DIST129:                RDIST=1.0D0/DIST
173:                DUMMY2=DFATT(RHO,DIST)*RDIST 
174:                CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)130:                CALL DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)
175:             ENDIF131:             ENDIF
176:          ENDDO132:          ENDDO
177:       ENDDO133:       ENDDO
178:    ENDDO134:    ENDDO
179: ENDDO135: ENDDO
180: 136: 
181: END SUBROUTINE DJWGR1137: END SUBROUTINE DJWGR1
182: 138: 
183: SUBROUTINE DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)139: SUBROUTINE DJWGR1GRAD(NATOMS,NPOS1,NPOS2,DUMMY2,RDIST,X,V)
184: USE MODHESS140: USE MODHESS
185: IMPLICIT NONE141: IMPLICIT NONE
186: INTEGER NPOS1, NPOS2, NATOMS, J1, J2, J3, J4142: INTEGER NPOS1, NPOS2, NATOMS, J1, J2, J3, J4
187: DOUBLE PRECISION X(3*NATOMS), DUMMY2, RDIST, DUMMY5, V(3*NATOMS)143: DOUBLE PRECISION X(3*NATOMS), DUMMY2, RDIST, DUMMY4, DUMMY5, V(3*NATOMS)
188: 144: 
 145: DUMMY4=DUMMY2*RDIST
189: J3=3*(NPOS1-1)146: J3=3*(NPOS1-1)
190: J4=3*(NPOS2-1)147: J4=3*(NPOS2-1)
191: 148: 
192: DO J1=1,3149: DO J1=1,3
193:    DUMMY5=DUMMY2*(X(J3+J1)-X(J4+J1))150:    DUMMY5=DUMMY4*(X(J3+J1)-X(J4+J1))
194:    V(J3+J1)=V(J3+J1)+DUMMY5151:    V(J3+J1)=V(J3+J1)+DUMMY5
195:    V(J4+J1)=V(J4+J1)-DUMMY5152:    V(J4+J1)=V(J4+J1)-DUMMY5
196: ENDDO153: ENDDO
197: 154: 
198: END SUBROUTINE DJWGR1GRAD155: END SUBROUTINE DJWGR1GRAD
199: 156: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0