hdiff output

r30227/djwgr1.f90 2016-04-04 22:30:04.236372898 +0100 r30226/djwgr1.f90 2016-04-04 22:30:04.624378003 +0100
 22: ! 22: !
 23: SUBROUTINE DJWGR1(NATOMS,X,V,ENERGY,GTEST,SECT) 23: SUBROUTINE DJWGR1(NATOMS,X,V,ENERGY,GTEST,SECT)
 24: USE MODHESS 24: USE MODHESS
 25: USE GENRIGID 25: USE GENRIGID
 26: IMPLICIT NONE 26: IMPLICIT NONE
 27: LOGICAL GTEST,SECT 27: LOGICAL GTEST,SECT
 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 FATT, DFATT, DDFATT, FREP, DFREP, DDFREP 30: DOUBLE PRECISION FATT, DFATT, DDFATT, FREP, DFREP, DDFREP
 31:  31: 
 32: ! 
 33: ! Derivatives of the pairwise site-site terms in terms of distance 
 34: ! 
 35: FATT(RHO,XDUMM)=-1.0D0 + (1.0D0 - EXP(RHO*(1.0D0 - XDUMM)))**2 32: 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 33: DFATT(RHO,XDUMM)=2.0D0*(-EXP(2.0D0*RHO*(1.0D0-XDUMM)) + EXP(RHO*(1.0D0-XDUMM)))*RHO
 37: DDFATT(RHO,XDUMM)=-2.0D0*(-2.0D0*EXP(2.0D0*RHO*(1.0D0-XDUMM)) + EXP(RHO*(1.0D0-XDUMM)))*RHO**2 34: DDFATT(RHO,XDUMM)=-2.0D0*(-2.0D0*EXP(2.0D0*RHO*(1.0D0-XDUMM)) + EXP(RHO*(1.0D0-XDUMM)))*RHO**2
 38: ! FATT(RHO,XDUMM)=0.0D0 
 39: ! DFATT(RHO,XDUMM)=0.0D0 
 40: ! DDFATT(RHO,XDUMM)=0.0D0 
 41: FREP(SIGMA,XDUMM)=(SIGMA/XDUMM)**12 35: FREP(SIGMA,XDUMM)=(SIGMA/XDUMM)**12
 42: DFREP(SIGMA,XDUMM)=-12.0D0*(SIGMA/XDUMM)**12/XDUMM 36: DFREP(SIGMA,XDUMM)=-12.0D0*(SIGMA/XDUMM)**12/XDUMM
 43: DDFREP(SIGMA,XDUMM)=156.0D0*(SIGMA/XDUMM)**12/XDUMM**2 37: DDFREP(SIGMA,XDUMM)=156.0D0*(SIGMA/XDUMM)**12/XDUMM**2
 44: ! FREP(SIGMA,XDUMM)=0.0D0 
 45: ! DFREP(SIGMA,XDUMM)=0.0D0 
 46: ! DDFREP(SIGMA,XDUMM)=0.0D0 
 47:  38: 
 48: ENERGY=0.0D0 39: ENERGY=0.0D0
 49: IF (GTEST) V(1:3*NATOMS)=0.0D0 40: IF (GTEST) V(1:3*NATOMS)=0.0D0
 50: IF (SECT) HESS(1:3*NATOMS,1:3*NATOMS)=0.0D0 41: IF (SECT) HESS(1:3*NATOMS,1:3*NATOMS)=0.0D0
 51:  42: 
 52: ! 
 53: ! 5 Morse plus two axial site pentamers from 
 54: ! S.N. Fejer, T. James, J. Hernandez-Rojas and D.J. Wales, Phys. Chem. Chem. Phys., 11, 2098-2104 (2009).  
 55: ! Energy Landscapes for Shells Assembled from Pentagonal and Hexagonal Pyramids  
 56: ! 
 57: RAD=5.0D0 43: RAD=5.0D0
 58: RHO=3.0D0 44: RHO=3.0D0
 59: SIGMA=(1.0D0+RAD*SQRT((5.0D0+SQRT(5.0D0))/2.0D0)) 45: SIGMA=(1.0D0+RAD*SQRT((5.0D0+SQRT(5.0D0))/2.0D0))
 60: EPSEFF=0.28D0 46: EPSEFF=0.28D0
 61:  47: 
 62: DO J1=1,NRIGIDBODY 48: DO J1=1,NRIGIDBODY
 63:    DO J2=J1+1,NRIGIDBODY 49:    DO J2=J1+1,NRIGIDBODY
 64: ! 
 65: ! Three different sorts of axial repulsion 
 66: ! 
 67:       NPOS1=RIGIDGROUPS(NSITEPERBODY(J1)-1,J1) 50:       NPOS1=RIGIDGROUPS(NSITEPERBODY(J1)-1,J1)
 68:       NPOS2=RIGIDGROUPS(NSITEPERBODY(J2)-1,J2) 51:       NPOS2=RIGIDGROUPS(NSITEPERBODY(J2)-1,J2)
 69:       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) 52:       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)
 70:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term 53:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term
 71:       IF (GTEST) THEN 54:       IF (GTEST) THEN
 72:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST) 55:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST)
 73:          RDIST=1.0D0/DIST 56:          RDIST=1.0D0/DIST
 74:          V(3*(NPOS1-1)+1)=V(3*(NPOS1-1)+1)+DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST 57:          V(3*(NPOS1-1)+1)=V(3*(NPOS1-1)+1)+DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST
 75:          V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST 58:          V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
 76:          V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST 59:          V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
 77:          V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST 60:          V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST
 78:          V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST 61:          V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
 79:          V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST 62:          V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
 80:       ENDIF 63:       ENDIF
 81:       IF (SECT) THEN 
 82:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST) 
 83:          DUMMY3=EPSEFF*DDFREP(SIGMA,DIST) 
 84:          RDIST=1.0D0/DIST 
 85: ! same rigid body, same Cartesian component 
 86:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)+ & 
 87:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST   
 88:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)+ & 
 89:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
 90:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)+ & 
 91:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
 92:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)+ & 
 93:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST   
 94:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)+ & 
 95:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
 96:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)+ & 
 97:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
 98: ! same rigid body, different Cartesian component 
 99:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)+ & 
100:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
101:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)+ & 
102:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
103:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)+ & 
104:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
105:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)+ & 
106:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
107:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)+ & 
108:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
109:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)+ & 
110:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
111:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)+ & 
112:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
113:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)+ & 
114:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
115:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)+ & 
116:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
117:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)+ & 
118:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
119:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)+ & 
120:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
121:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)+ & 
122:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
123: ! different rigid body, same Cartesian component 
124:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)- & 
125:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
126:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)- & 
127:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
128:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)- & 
129:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
130:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)- & 
131:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
132:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)- & 
133:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
134:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)- & 
135:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
136: ! different rigid body, different Cartesian component 
137:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)- & 
138:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
139:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)- & 
140:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
141:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)- & 
142:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
143:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)- & 
144:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
145:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)- & 
146:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
147:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)- & 
148:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
149:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)- & 
150:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
151:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)- & 
152:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
153:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)- & 
154:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
155:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)- & 
156:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
157:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)- & 
158:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
159:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)- & 
160:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
161:       ENDIF 
162:       NPOS1=RIGIDGROUPS(NSITEPERBODY(J1)-1,J1) 64:       NPOS1=RIGIDGROUPS(NSITEPERBODY(J1)-1,J1)
163:       NPOS2=RIGIDGROUPS(NSITEPERBODY(J2),J2) 65:       NPOS2=RIGIDGROUPS(NSITEPERBODY(J2),J2)
164:       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) 66:       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)
165:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term 67:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial-axial repulsive term
166:       IF (GTEST) THEN 68:       IF (GTEST) THEN
167:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST) 69:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST)
168:          RDIST=1.0D0/DIST 70:          RDIST=1.0D0/DIST
169:          V(3*(NPOS1-1)+1)=V(3*(NPOS1-1)+1)+DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST 71:          V(3*(NPOS1-1)+1)=V(3*(NPOS1-1)+1)+DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST
170:          V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST 72:          V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
171:          V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST 73:          V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
172:          V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST 74:          V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST
173:          V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST 75:          V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
174:          V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST 76:          V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
175:       ENDIF 77:       ENDIF
176:       IF (SECT) THEN 
177:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST) 
178:          DUMMY3=EPSEFF*DDFREP(SIGMA,DIST) 
179:          RDIST=1.0D0/DIST 
180: ! same rigid body, same Cartesian component 
181:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)+ & 
182:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST   
183:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)+ & 
184:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
185:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)+ & 
186:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
187:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)+ & 
188:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST   
189:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)+ & 
190:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
191:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)+ & 
192:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST 
193: ! same rigid body, different Cartesian component 
194:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)+ & 
195:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
196:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)+ & 
197:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
198:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)+ & 
199:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
200:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)+ & 
201:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
202:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)+ & 
203:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
204:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)+ & 
205:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
206:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)+ & 
207:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
208:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)+ & 
209:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
210:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)+ & 
211:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
212:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)+ & 
213:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
214:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)+ & 
215:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
216:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)+ & 
217:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
218: ! different rigid body, same Cartesian component 
219:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)- & 
220:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
221:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)- & 
222:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
223:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)- & 
224:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
225:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)- & 
226:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
227:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)- & 
228:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
229:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)- & 
230:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
231: ! different rigid body, different Cartesian component 
232:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)- & 
233:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
234:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)- & 
235:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
236:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)- & 
237:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
238:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)- & 
239:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
240:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)- & 
241:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
242:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)- & 
243:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
244:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)- & 
245:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
246:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)- & 
247:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
248:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)- & 
249:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
250:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)- & 
251:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
252:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)- & 
253:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
254:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)- & 
255:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
256:       ENDIF 
257:       NPOS1=RIGIDGROUPS(NSITEPERBODY(J1),J1) 78:       NPOS1=RIGIDGROUPS(NSITEPERBODY(J1),J1)
258:       NPOS2=RIGIDGROUPS(NSITEPERBODY(J2)-1,J2) 79:       NPOS2=RIGIDGROUPS(NSITEPERBODY(J2)-1,J2)
259:       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) 80:       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)
260:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial 2-axial repulsive term 81:       ENERGY=ENERGY+EPSEFF*FREP(SIGMA,DIST)  ! axial 2-axial repulsive term
261:       IF (GTEST) THEN 82:       IF (GTEST) THEN
262:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST) 83:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST)
263:          RDIST=1.0D0/DIST 84:          RDIST=1.0D0/DIST
264:          V(3*(NPOS1-1)+1)=V(3*(NPOS1-1)+1)+DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST 85:          V(3*(NPOS1-1)+1)=V(3*(NPOS1-1)+1)+DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST
265:          V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST 86:          V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
266:          V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST 87:          V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
267:          V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST 88:          V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST
268:          V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST 89:          V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
269:          V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST 90:          V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
270:       ENDIF 91:       ENDIF
271:       IF (SECT) THEN 92:       IF (SECT) THEN
272:          DUMMY2=EPSEFF*DFREP(SIGMA,DIST) 93:          DUMMY2=DFATT(RHO,DIST)
273:          DUMMY3=EPSEFF*DDFREP(SIGMA,DIST) 94:          DUMMY3=DDFATT(RHO,DIST)
274:          RDIST=1.0D0/DIST 95:          RDIST=1.0D0/DIST
275: ! same rigid body, same Cartesian component 96: !        HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=
276:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)+ & 97: !        HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=
277:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST   98: !        HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=
278:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)+ & 99: !        HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=
279:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST100: !        HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=
280:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)+ &101: !        HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=
281:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST102: !        HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=
282:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)+ &103: !        HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=
283:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST  104: !        HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=
284:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)+ &105: !        HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=
285:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST106: !        HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=
286:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)+ &107: !        HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=
287:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST108: !        HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=
288: ! same rigid body, different Cartesian component109: !        HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=
289:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)+ &110: !        HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=
290:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)   111: !        HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=
291:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)+ &112: !        HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=
292:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)113: !        HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=
293:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)+ & 
294:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
295:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)+ & 
296:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
297:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)+ & 
298:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
299:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)+ & 
300:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
301:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)+ & 
302:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
303:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)+ & 
304:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
305:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)+ & 
306:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
307:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)+ & 
308:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
309:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)+ & 
310:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
311:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)+ & 
312:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
313: ! different rigid body, same Cartesian component 
314:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)- & 
315:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
316:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)- & 
317:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
318:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)- & 
319:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
320:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)- & 
321:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
322:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)- & 
323:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
324:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)- & 
325:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
326: ! different rigid body, different Cartesian component 
327:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)- & 
328:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
329:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)- & 
330:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
331:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)- & 
332:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
333:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)- & 
334:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
335:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)- & 
336:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
337:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)- & 
338:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
339:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)- & 
340:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
341:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)- & 
342:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
343:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)- & 
344:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
345:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)- & 
346:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
347:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)- & 
348:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
349:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)- & 
350:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
351:       ENDIF114:       ENDIF
352: !115: !
353: ! Sum over the attractive sites116: ! Sum over the attractive sites
354: !117: !
355:       DO J3=1,NSITEPERBODY(J1)-2     ! # Morse sites in rb J1 NSITEPERBODY(J1)118:       DO J3=1,NSITEPERBODY(J1)-2     ! # Morse sites in rb J1 NSITEPERBODY(J1)
356:          NPOS1=RIGIDGROUPS(J3,J1)    ! where is this site in the list?119:          NPOS1=RIGIDGROUPS(J3,J1)    ! where is this site in the list?
357:          DO J4=1,NSITEPERBODY(J2)-2  ! # Morse sites in rb J2 NSITEPERBODY(J2)120:          DO J4=1,NSITEPERBODY(J2)-2  ! # Morse sites in rb J2 NSITEPERBODY(J2)
358:             NPOS2=RIGIDGROUPS(J4,J2) ! where is this site in the list?121:             NPOS2=RIGIDGROUPS(J4,J2) ! where is this site in the list?
359:             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)122:             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)
360:             ENERGY=ENERGY+FATT(RHO,DIST) 123:             ENERGY=ENERGY+FATT(RHO,DIST) 
365:                V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST128:                V(3*(NPOS1-1)+2)=V(3*(NPOS1-1)+2)+DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
366:                V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST129:                V(3*(NPOS1-1)+3)=V(3*(NPOS1-1)+3)+DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
367:                V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST130:                V(3*(NPOS2-1)+1)=V(3*(NPOS2-1)+1)-DUMMY2*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST
368:                V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST131:                V(3*(NPOS2-1)+2)=V(3*(NPOS2-1)+2)-DUMMY2*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST
369:                V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST132:                V(3*(NPOS2-1)+3)=V(3*(NPOS2-1)+3)-DUMMY2*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST
370:             ENDIF133:             ENDIF
371:             IF (SECT) THEN134:             IF (SECT) THEN
372:                DUMMY2=DFATT(RHO,DIST)135:                DUMMY2=DFATT(RHO,DIST)
373:                DUMMY3=DDFATT(RHO,DIST)136:                DUMMY3=DDFATT(RHO,DIST)
374:                RDIST=1.0D0/DIST137:                RDIST=1.0D0/DIST
375: ! same rigid body, same Cartesian component138: !              HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=
376:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+1)+ &139: !              HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=
377:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST  140: !              HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=
378:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+2)+ &141: !              HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=
379:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST142: !              HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=
380:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+3)+ &143: !              HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=
381:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST144: !              HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=
382:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+1)+ &145: !              HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=
383:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST  146: !              HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=
384:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+2)+ &147: !              HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=
385:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST148: !              HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=
386:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+3)+ &149: !              HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=
387:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)+DUMMY2*RDIST150: !              HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=
388: ! same rigid body, different Cartesian component151: !              HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=
389:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+2)+ &152: !              HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=
390:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)   153: !              HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=
391:          HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS1-1)+3)+ &154: !              HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=
392:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)155: !              HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=
393:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+1)+ & 
394:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
395:          HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS1-1)+3)+ & 
396:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
397:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+1)+ & 
398:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
399:          HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS1-1)+2)+ & 
400:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
401:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+2)+ & 
402:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
403:          HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS2-1)+3)+ & 
404:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
405:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+1)+ & 
406:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
407:          HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS2-1)+3)+ & 
408:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
409:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+1)+ & 
410:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
411:          HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS2-1)+2)+ & 
412:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
413: ! different rigid body, same Cartesian component 
414:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+1)- & 
415:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
416:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+2)- & 
417:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
418:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+3)- & 
419:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
420:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+1)- & 
421:   &                                  ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
422:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+2)- & 
423:   &                                  ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
424:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+3)- & 
425:   &                                  ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*RDIST)**2*(DUMMY3-DUMMY2*RDIST)-DUMMY2*RDIST 
426: ! different rigid body, different Cartesian component 
427:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+2)- & 
428:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
429:          HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+1,3*(NPOS2-1)+3)- & 
430:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
431:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+1)- & 
432:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
433:          HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)=HESS(3*(NPOS1-1)+2,3*(NPOS2-1)+3)- & 
434:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
435:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+1)- & 
436:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
437:          HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)=HESS(3*(NPOS1-1)+3,3*(NPOS2-1)+2)- & 
438:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
439:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+2)- & 
440:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST)    
441:          HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+1,3*(NPOS1-1)+3)- & 
442:   &                         ((X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
443:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+1)- & 
444:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
445:          HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)=HESS(3*(NPOS2-1)+2,3*(NPOS1-1)+3)- & 
446:   &                         ((X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2))*(X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
447:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+1)- & 
448:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+1)-X(3*(NPOS2-1)+1)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
449:          HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)=HESS(3*(NPOS2-1)+3,3*(NPOS1-1)+2)- & 
450:   &                         ((X(3*(NPOS1-1)+3)-X(3*(NPOS2-1)+3))*(X(3*(NPOS1-1)+2)-X(3*(NPOS2-1)+2)))*RDIST**2*(DUMMY3-DUMMY2*RDIST) 
451:             ENDIF156:             ENDIF
452:          ENDDO157:          ENDDO
453:       ENDDO158:       ENDDO
454:    ENDDO159:    ENDDO
455: ENDDO160: ENDDO
456: 161: 
457: END SUBROUTINE DJWGR1162: END SUBROUTINE DJWGR1


r30227/potential.f 2016-04-04 22:30:04.432375479 +0100 r30226/potential.f 2016-04-04 22:30:04.820380580 +0100
1630:             CALL MULTISITEPY2 (COORDS, VNEW, ENERGY, GTEST)1630:             CALL MULTISITEPY2 (COORDS, VNEW, ENERGY, GTEST)
1631:             IF (SSTEST) THEN1631:             IF (SSTEST) THEN
1632:                CALL MULTISITEPYSECDER(COORDS,SSTEST)1632:                CALL MULTISITEPYSECDER(COORDS,SSTEST)
1633:             END IF1633:             END IF
1634:             IF (PTEST) THEN1634:             IF (PTEST) THEN
1635:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1635:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1636:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1636:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1637:             END IF1637:             END IF
1638:          ELSE IF (DJWRBT) THEN1638:          ELSE IF (DJWRBT) THEN
1639:             IF (DJWRBID.EQ.1) THEN1639:             IF (DJWRBID.EQ.1) THEN
1640: !              CALL DJWGR1(NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST)1640:                CALL DJWGR1(NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST)
1641:  
1642:                DIFF=1.0D-4 
1643:                PRINT*,'analytic and numerical gradients:' 
1644:                IF (.NOT.(ALLOCATED(HESS))) ALLOCATE(HESS(3*NATOMS,3*NATOMS)) 
1645:                CALL DJWGR1(NATOMS,COORDS, VNEW, ENERGY, .TRUE., .TRUE.) 
1646:                PRINT '(A,I8)','SIZE(HESS)=',SIZE(HESS) 
1647:                HESSDUM(1:3*NATOMS,1:3*NATOMS)=HESS(1:3*NATOMS,1:3*NATOMS) 
1648:                DO J1=1,3*NATOMS 
1649:                   COORDS(J1)=COORDS(J1)+DIFF 
1650:                   CALL DJWGR1(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.) 
1651:                   COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
1652:                   CALL DJWGR1(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.) 
1653:                   COORDS(J1)=COORDS(J1)+DIFF 
1654:                   WRITE(*,'(I5,2F20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF) 
1655:                ENDDO 
1656:                PRINT*,'analytic and numerical second derivatives:' 
1657:                DO J1=1,3*NATOMS 
1658:                   COORDS(J1)=COORDS(J1)+DIFF 
1659:                   CALL DJWGR1(NATOMS,COORDS,VPLUS,EPLUS,.TRUE.,.FALSE.) 
1660:                   COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
1661:                   CALL DJWGR1(NATOMS,COORDS,VMINUS,EMINUS,.TRUE.,.FALSE.) 
1662:                   COORDS(J1)=COORDS(J1)+DIFF 
1663:                   DO J2=1,3*NATOMS 
1664:                      DUMMY1=(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF) 
1665:                      IF ((ABS(DUMMY1).GT.1.0D-10).AND.  
1666:      &                   (ABS(100.0D0*(HESS(J1,J2)-DUMMY1)/DUMMY1).GT.1.0D0)) THEN 
1667:                      WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),DUMMY1,'   X' 
1668:                      ELSE 
1669:                         WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),DUMMY1 
1670:                      ENDIF 
1671:                   ENDDO 
1672:                ENDDO 
1673:                STOP 
1674:             ENDIF1641:             ENDIF
1675:             IF (PTEST) THEN1642:             IF (PTEST) THEN
1676:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1643:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1677:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1644:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1678:             ENDIF1645:             ENDIF
1679:          ELSE IF (NCAPT) THEN1646:          ELSE IF (NCAPT) THEN
1680:             CALL NEWCAPSID (COORDS,VNEW,ENERGY,GTEST,SSTEST)1647:             CALL NEWCAPSID (COORDS,VNEW,ENERGY,GTEST,SSTEST)
1681:             IF (PTEST) THEN1648:             IF (PTEST) THEN
1682:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1649:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1683:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1650:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0