hdiff output

r30590/SBM.f 2016-06-14 00:30:08.359777601 +0100 r30589/SBM.f 2016-06-14 00:30:08.603780883 +0100
  4:       INTEGER NATOMS,i,j  4:       INTEGER NATOMS,i,j
  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)
  6:       DOUBLE PRECISION ENERGY,STT  6:       DOUBLE PRECISION ENERGY,STT
  7:   7: 
  8:       LOGICAL :: CALLED=.FALSE.  8:       LOGICAL :: CALLED=.FALSE.
  9:       LOGICAL GTEST, STEST  9:       LOGICAL GTEST, STEST
 10:         integer NSBMMAX 10:         integer NSBMMAX
 11:         parameter(NSBMMAX=5000) 11:         parameter(NSBMMAX=5000)
 12:       DOUBLE PRECISION  Rb(NSBMMAX), bK(NSBMMAX), ANTC(NSBMMAX), 12:       DOUBLE PRECISION  Rb(NSBMMAX), bK(NSBMMAX), ANTC(NSBMMAX),
 13:      Q Tk(NSBMMAX), PK(NSBMMAX), PHISBM(NSBMMAX), Sigma(NSBMMAX*5), 13:      Q Tk(NSBMMAX), PK(NSBMMAX), PHISBM(NSBMMAX), Sigma(NSBMMAX*5),
 14:      Q EpsC(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 14:      Q EpsC(NSBMMAX*5),NNCsigma(NSBMMAX,NSBMMAX),NCswitch,NCcut
 15:       INTEGER  Ib1(NSBMMAX), Ib2(NSBMMAX), IT(NSBMMAX), JT(NSBMMAX), 15:       INTEGER  Ib1(NSBMMAX), Ib2(NSBMMAX), IT(NSBMMAX), JT(NSBMMAX),
 16:      Q KT(NSBMMAX),IP(NSBMMAX), JP(NSBMMAX), KP(NSBMMAX), 16:      Q KT(NSBMMAX),IP(NSBMMAX), JP(NSBMMAX), KP(NSBMMAX),
 17:      Q LP(NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 17:      Q LP(NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5),
 18:      Q PHITYPE(NSBMMAX),NBA,NTA,NPA,NC,NNCinc(NSBMMAX,NSBMMAX),SBMPRN, SBMPRI(NSBMMAX) 18:      Q  PHITYPE(NSBMMAX),NBA, NTA, NPA, NC
 19:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3) 19:         integer CONTACTTYPE
 20:         integer CONTACTTYPE,SBMCHARGEON(NSBMMAX) 20:         logical used(NSBMMAX,NSBMMAX)
 21:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 21:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 22:      Q sigma, epsC,NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 22:      Q sigma, epsC, NNCsigma,NCswitch,NCcut,STT
 23:      Q SBMPRK,SBMPRX 
 24:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP, 23:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP,
 25:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC,NNCinc, 24:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC,
 26:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI 25:      Q  CONTACTTYPE
  26:         common /logical/ used
  27: 
 27:         if(NATOMS.gt. NSBMMAX)then 28:         if(NATOMS.gt. NSBMMAX)then
 28:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 29:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 29:         STOP 30:         STOP
 30:         endif 31:         endif
 31:  32: 
 32:        if(.NOT.CALLED)then 33:        if(.NOT.CALLED)then
 33:         write(*,*) 34: 
 34:         write(*,*)  'Calculations will use a Structure-based SMOG model, described in:' 35:         write(*,*)  'Using a Structure-based SMOG model, described in:'
 35:         write(*,*)  'Whitford, et al. Prot. Struct. Func. Bioinfo. 75, 430-441, 2009.' 36:         write(*,*)  'Whitford, et al. Prot. Struct. Func. Bioinfo. 75, 430-441, 2009.'
 36:         write(*,*) 37: 
 37:          38:           do i=1,NATOMS-1
  39:           do j=i+1,NATOMS
  40:                 used(i,j) =.FALSE.
  41:                 used(j,i) =.FALSE.
  42:           enddo
  43:           enddo
  44: 
  45: 
  46: 
 38:  47: 
 39:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP,  48:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP, 
 40:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, Sigma,  49:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, Sigma, 
 41:      Q EpsC, NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT, 50:      Q EpsC, NNCsigma,NBA, NTA, NPA, NC, NCswitch,NCcut,STT,
 42:      Q NSBMMAX,CONTACTTYPE, 51:      Q NSBMMAX,CONTACTTYPE)
 43:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 52: 
 44:         if(NBA .gt. NSBMMAX .or. NTA .gt. NSBMMAX .or. 53:         if(NBA .gt. NSBMMAX .or. NTA .gt. NSBMMAX .or.
 45:      Q  NPA .gt. NSBMMAX .or. NC .gt. 5*NSBMMAX)then 54:      Q  NPA .gt. NSBMMAX .or. NC .gt. 5*NSBMMAX)then
 46:         write(*,*) 'increase array size' 55:         write(*,*) 'increase array size'
 47:         stop 56:         stop
 48:         endif  57:         endif 
 49:  58: 
 50:         CALLED=.TRUE. 59:         CALLED=.TRUE.
 51:         endIF 60:         endIF
 52: ! call the energy routine 61: ! call the energy routine
 53:  62: 
 54:       call calc_energy_SBM(qo,natoms, GRAD, energy, Ib1, Ib2, 63:       call calc_energy_SBM(qo,natoms, GRAD, energy, Ib1, Ib2,
 55:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK, 64:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,
 56:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC,  65:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC, 
 57:      Q  NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON, 66:      Q  NNCsigma,NBA, NTA, NPA, NC, NCswitch,STT,NCcut,CONTACTTYPE,used)
 58:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 67: 
 59:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 
 60:  68: 
 61:       IF (STEST) THEN 69:       IF (STEST) THEN
 62:          PRINT '(A)','ERROR - second derivatives not available' 70:          PRINT '(A)','ERROR - second derivatives not available'
 63:          STOP 71:          STOP
 64:       ENDIF 72:       ENDIF
 65:       return 73:       return
 66:       end 74:       end
 67:  75: 
 68:  76: 
 69:  77: 
 70:  78: 
 71:  79: 
 72: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 80: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 73: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 81: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 74: !* startt then the velocities are assigned, otherwise, they are read   * 82: !* startt then the velocities are assigned, otherwise, they are read   *
 75: !* by selecting 2, or generated by selecting 3                         * 83: !* by selecting 2, or generated by selecting 3                         *
 76: !*********************************************************************** 84: !***********************************************************************
 77:  85: 
 78:       subroutine SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, 86:       subroutine SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk,
 79:      Q  IP, JP, KP, LP, PK,PHITYPE,PHISBM,IC, JC, Sigma,  87:      Q  IP, JP, KP, LP, PK,PHITYPE,PHISBM,IC, JC, Sigma, 
 80:      Q EpsC, NNCinc, NBA,  88:      Q EpsC, NNCsigma, NBA, 
 81:      Q NTA, NPA, NC, SBMCHARGE, SBMCHARGEON,NCswitch,NCcut,STT,NSBMMAX, 89:      Q NTA, NPA, NC,NCswitch,NCcut,STT,NSBMMAX,CONTACTTYPE)
 82:      Q CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 
 83:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 
 84:       USE KEY 90:       USE KEY
 85:       USE COMMONS, only: ATMASS 
 86:       implicit NONE 91:       implicit NONE
 87:  92: 
 88:         integer i,j,MaxCon,NNCmax,NATOMS,storage, dummy,  ANr, IB11, 93:         integer i,j,MaxCon,NNCmax,NATOMS,storage, dummy,  ANr, IB11,
 89:      Q IB12, Ib22, Ib21,IT1, JT1, KT1, IT2, JT2, KT2, IP1, JP1, 94:      Q IB12, Ib22, Ib21,IT1, JT1, KT1, IT2, JT2, KT2, IP1, JP1,
 90:      Q KP1, LP1, IP2, JP2, KP2, 95:      Q KP1, LP1, IP2, JP2, KP2,
 91:      Q LP2, nBA1, nTA1, nPA1, nBA2, nTA2, nPA2,  ind1, ind2, ANt, 96:      Q LP2, nBA1, nTA1, nPA1, nBA2, nTA2, nPA2,  ind1, ind2, ANt,
 92:      Q  MDT1, MDT2, cl1, cl2,tempi,NSBMMAX 97:      Q  MDT1, MDT2, cl1, cl2,tempi,NSBMMAX
 93:  98: 
 94:       DOUBLE PRECISION  Rb(NSBMMAX), bK(NSBMMAX), ANTC(NSBMMAX),  99:       DOUBLE PRECISION  Rb(NSBMMAX), bK(NSBMMAX), ANTC(NSBMMAX), 
 95:      Q Tk(NSBMMAX), PK(NSBMMAX), PHISBM(NSBMMAX),100:      Q Tk(NSBMMAX), PK(NSBMMAX), PHISBM(NSBMMAX),
 96:      Q Sigma(NSBMMAX*5), EpsC(NSBMMAX*5),  101:      Q Sigma(NSBMMAX*5), EpsC(NSBMMAX*5),  
 97:      Q NCswitch,NCcut,STT,SBMCHARGE(NATOMS)102:      Q NNCsigma(NATOMS,NATOMS),NCswitch,NCcut,STT
 98:       INTEGER  Ib1(NSBMMAX), Ib2(NSBMMAX), IT(NSBMMAX), JT(NSBMMAX), 103:       INTEGER  Ib1(NSBMMAX), Ib2(NSBMMAX), IT(NSBMMAX), JT(NSBMMAX), 
 99:      Q KT(NSBMMAX),IP(NSBMMAX), JP(NSBMMAX), KP(NSBMMAX),104:      Q KT(NSBMMAX),IP(NSBMMAX), JP(NSBMMAX), KP(NSBMMAX),
100:      Q LP(NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 105:      Q LP(NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 
101:      Q  NBA, NTA, NPA, NC,SBMCHARGEON(NATOMS)106:      Q  NBA, NTA, NPA, NC
102:       INTEGER PHITYPE(NSBMMAX),NNCinc(NATOMS,NATOMS)107:       INTEGER PHITYPE(NSBMMAX)
103:        DOUBLE PRECISION  pinitmax, TK1, TK2,  APTtemp, msT,108:        DOUBLE PRECISION  pinitmax, TK1, TK2,  APTtemp, msT,
104:      Q SigmaT1, SigmaT2, epstemp109:      Q SigmaT1, SigmaT2, epstemp
 110:         character(LEN=20) FMTB, FMTT, FMTP, CA, RP
105:       integer AA,BB,ANTEMP111:       integer AA,BB,ANTEMP
106:       DOUBLE PRECISION dx,dy,dz112:       DOUBLE PRECISION dx,dy,dz
107:        double precision PI113:        double precision PI
108:         DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut114:         DOUBLE PRECISION RSig, Reps
109:         logical TEMPARRAY115:         logical TEMPARRAY
110:         dimension TEMPARRAY(NATOMS,NATOMS)116:         dimension TEMPARRAY(NATOMS,NATOMS)
111:         integer CONTACTTYPE117:         integer CONTACTTYPE
112:         character TMPCHAR 
113:         integer TMPINT,NUMCHARGES,nexc,I1,I2 
114:         double precision TMPREAL,concentration 
115:  
116:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS) 
117:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3) 
118:  
119:  
120:  
121:       pi = 3.14159265358979323846264338327950288419716939937510118:       pi = 3.14159265358979323846264338327950288419716939937510
122: 119: 
123:         NNCmax = NATOMS*NATOMS120:         NNCmax = NATOMS*NATOMS
124:         MaxCon=NATOMS*5121:         MaxCon=NATOMS*5
 122: ! old formatting
 123:         FMTB="(3I5,2F8.3)"
 124:         FMTT="(4I5,2F8.3)"
 125:         FMTP="(5I5,2F8.3)"
 126:         CA="(3I5,F10.3, F9.6)"
 127:         RP="(I5,2I5, 2F8.3)"
125: 128: 
126:         do i=1,NATOMS129:         do i=1,NATOMS
127:                 do j=1,NATOMS130:                 do j=1,NATOMS
128:                 TEMPARRAY(i,j)=.TRUE.131:                 TEMPARRAY(i,j)=.TRUE.
129:                 enddo132:                 enddo
130:         enddo133:         enddo
131: 134: 
132: 135: 
133: 136: 
134: ! These lines read in the parameters.137: ! These lines read in the parameters.
135:         open(30, file='SBM.INP', status='old', access='sequential')138:         open(30, file='SBM.INP', status='old', access='sequential')
136:         write(*,*) 139: 
137:         write(*,*) 'Reading the forcefield from SBM.INP' 
138:         write(*,*)  
139:         read(30,*) 
140:         read(30,*)140:         read(30,*)
141:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut 
142:         read(30,*)141:         read(30,*)
143:         read(30,*) RSig, Reps,NCswitch,NCcut142:         read(30,*) RSig, Reps,NCswitch,NCcut
144:         write(*,*) 'RSig, Reps,NCswitch,NCcut' 
145:         write(*,*) RSig, Reps,NCswitch,NCcut143:         write(*,*) RSig, Reps,NCswitch,NCcut
146:         write(*,*) 'Reading electrostatic parameters'144:         
147:         write(*,*) PREFACTOR,DC,KAPPA,DHswitch,DHcut 
148:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is is units 
149:         ! A^-1 
150:         KAPPA=0.32*sqrt(CONCENTRATION) 
151:         PREFACTOR=PREFACTOR/DC 
152:         read(30,*) ANtemp145:         read(30,*) ANtemp
153:         write(*,*) ANtemp, ' atoms'146:         write(*,*) ANtemp
154:         NUMCHARGES=1 
155:         do i=1, ANtemp147:         do i=1, ANtemp
156:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)148:           read(30,*) 
157:           if(SBMCHARGE(i) .ne. 0)then 
158:              NUMCHARGES=NUMCHARGES+1 
159:              SBMCHARGEON(NUMCHARGES)=i 
160:           endif  
161:         end do149:         end do
162:         SBMCHARGEON(1)=NUMCHARGES-1150: 
163:          
164:         read(30,*) 151:         read(30,*) 
165:         read(30,*) NC152:         read(30,*) NC
166:         write(*,*) NC, 'contacts'        153:         write(*,*) NC, 'contacts'        
167:         read(30,*) CONTACTTYPE154:         read(30,*) CONTACTTYPE
 155: 
168:           if(NC .gt. MaxCon)then156:           if(NC .gt. MaxCon)then
169:              write(*,*) 'too many contacts'157:              write(*,*) 'too many contacts'
170:              STOP158:              STOP
171:           endif159:           endif
172: 160: 
173:         do i=1, NC161:         do i=1, NC
174:           read(30, *) IC(i), JC(i), Sigma(i), EpsC(i)162:           read(30, *) IC(i), JC(i), Sigma(i), EpsC(i)
175:           TEMPARRAY(IC(i), JC(i))=.FALSE. 
176:           TEMPARRAY(JC(i), IC(i))=.FALSE. 
177:           Sigma(i)=Sigma(i)**2 
178:         end do163:         end do
179: 164: 
180:           read(30,*)165:           read(30,*)
181:           read(30,*) nBA166:           read(30,*) nBA
182:           write(*,*) nBA, 'bonds'167: 
183:         do i=1, nBA168:         do i=1, nBA
184:           read(30,*) Ib1(i), Ib2(i),Rb(i), bK(i)169:           read(30,*) Ib1(i), Ib2(i),Rb(i), bK(i)
185:           TEMPARRAY(Ib1(i), Ib2(i))=.FALSE.170:           TEMPARRAY(Ib1(i), Ib2(i))=.FALSE.
186:           TEMPARRAY(Ib2(i), Ib1(i))=.FALSE.171:           TEMPARRAY(Ib2(i), Ib1(i))=.FALSE.
187:         end do172:         end do
188: 173: 
 174: ! read non-native interactions
 175: !        read(30,*)
 176: !        read(30,*) NNC
 177: 
 178: !        do i=1,NATOMS
 179: !                do j=1,NATOMS
 180: !                TEMPARRAY(i,j)=.TRUE.
 181: !                enddo
 182: !        enddo
 183: 
 184: !        if(NNC .gt. NNCmax)then
 185: !        write(*,*) 'too many non contacts'
 186: !        STOP
 187: !        endif        
 188: !        do i=1, NNC
 189: !           read(30,*) AA,BB
 190: !           TEMPARRAY(AA,BB)=.FALSE.
 191: !           TEMPARRAY(BB,AA)=.FALSE.
 192: !        enddo
 193: 
 194: !        tempi=0
 195: 
 196: !        do i=1,NATOMS-1
 197: !                do j=i+1,NATOMS
 198: !                if(TEMPARRAY(i,j) .eqv. .TRUE.)then
 199: !                  tempi=tempi+1
 200: !                  INC(tempi)=i
 201: !                  JNC(tempi)=j
 202: !                  NCsigma(tempi)=rsig**2
 203: !                  NNCeps(tempi)=reps
 204: ! this simplifies calculations later
 205: !         NNCsigma(tempi) = NNCEps(tempi)*NCsigma(tempi)**6
 206: !                endif
 207: !                enddo
 208: !        end do
 209: 
 210: !        NNC=tempi
 211: !        write(*,*) NNC, 'noc'
 212: 
189:           read(30,*)213:           read(30,*)
190:           read(30,*) nTA214:           read(30,*) nTA
191:           write(*,*) nTA, 'angles' 
192:         do i=1, nTA215:         do i=1, nTA
193:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)216:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)
194:           TEMPARRAY(IT(i), JT(i))=.FALSE.217:           TEMPARRAY(IT(i), JT(i))=.FALSE.
195:           TEMPARRAY(JT(i), IT(i))=.FALSE.218:           TEMPARRAY(JT(i), IT(i))=.FALSE.
196:           TEMPARRAY(IT(i), KT(i))=.FALSE.219:           TEMPARRAY(IT(i), KT(i))=.FALSE.
197:           TEMPARRAY(KT(i), IT(i))=.FALSE.220:           TEMPARRAY(KT(i), IT(i))=.FALSE.
198:           TEMPARRAY(JT(i), KT(i))=.FALSE.221:           TEMPARRAY(JT(i), KT(i))=.FALSE.
199:           TEMPARRAY(KT(i), JT(i))=.FALSE.222:           TEMPARRAY(KT(i), JT(i))=.FALSE.
200:         enddo223:         enddo
201: 224: 
202:           read(30,*) 225:           read(30,*) 
203:           read(30,*) nPA226:           read(30,*) nPA
204:           write(*,*) nPA, 'dihedrals' 
205: 227: 
206: ! this reads in the dihedral angles and calculates the cosines and sines228: ! this reads in the dihedral angles and calculates the cosines and sines
207: ! in order to make the force and energy calculations easier, later.229: ! in order to make the force and energy calculations easier, later.
208:         do i=1, npA230:         do i=1, npA
209:            read(30,*) IP(i),JP(i),KP(i),LP(i),231:            read(30,*) IP(i),JP(i),KP(i),LP(i),
210:      Q  PHITYPE(i),PHISBM(i),PK(i)232:      Q  PHITYPE(i),PHISBM(i),PK(i)
211:           TEMPARRAY(IP(i), JP(i))=.FALSE.233:           TEMPARRAY(IP(i), JP(i))=.FALSE.
212:           TEMPARRAY(JP(i), IP(i))=.FALSE.234:           TEMPARRAY(JP(i), IP(i))=.FALSE.
213:           TEMPARRAY(IP(i), KP(i))=.FALSE.235:           TEMPARRAY(IP(i), KP(i))=.FALSE.
214:           TEMPARRAY(KP(i), IP(i))=.FALSE.236:           TEMPARRAY(KP(i), IP(i))=.FALSE.
216:           TEMPARRAY(LP(i), IP(i))=.FALSE.238:           TEMPARRAY(LP(i), IP(i))=.FALSE.
217:           TEMPARRAY(JP(i), KP(i))=.FALSE.239:           TEMPARRAY(JP(i), KP(i))=.FALSE.
218:           TEMPARRAY(KP(i), JP(i))=.FALSE.240:           TEMPARRAY(KP(i), JP(i))=.FALSE.
219:           TEMPARRAY(JP(i), LP(i))=.FALSE.241:           TEMPARRAY(JP(i), LP(i))=.FALSE.
220:           TEMPARRAY(LP(i), JP(i))=.FALSE.242:           TEMPARRAY(LP(i), JP(i))=.FALSE.
221:           TEMPARRAY(KP(i), LP(i))=.FALSE.243:           TEMPARRAY(KP(i), LP(i))=.FALSE.
222:           TEMPARRAY(LP(i), KP(i))=.FALSE.244:           TEMPARRAY(LP(i), KP(i))=.FALSE.
223:         245:         
224:         END DO246:         END DO
225: 247: 
226:           read(30,*) 
227:           read(30,*) nexc 
228:           write(*,*) nexc, 'exculusions' 
229:  
230:         do i=1, nexc 
231:           read(30,*) I1, I2 
232:           TEMPARRAY(I1, I2)=.FALSE. 
233:           TEMPARRAY(I2, I1)=.FALSE. 
234:         end do 
235:  
236: 248: 
237:         do i=1,NATOMS-1249:         do i=1,NATOMS-1
238:                 do j=i+1,NATOMS250:                 do j=i+1,NATOMS
239:                 if(TEMPARRAY(i,j) .eqv. .TRUE.)then251:                 if(TEMPARRAY(i,j) .eqv. .TRUE.)then
240: ! this simplifies calculations later252: ! this simplifies calculations later
241:          NNCinc(i,j) = 1253:          NNCsigma(i,j) = reps*rsig**12
242:          NNCinc(j,i) = 0254:          NNCsigma(j,i) = reps*rsig**12
243:                 else255:                 else
244:          NNCinc(i,j) = 0256:          NNCsigma(i,j) = 0
245:          NNCinc(j,i) = 0257:          NNCsigma(j,i) = 0
246: 258: 
247:                 endif259:                 endif
248:                 enddo260:                 enddo
249:         end do261:         end do
250: 262: 
251:                 STT=reps*rsig**12263:                 STT=reps*rsig**12
252: 264: !        read(30,*) AN
253:         do I=1,NATOMS 
254:           TARR(I)=0 
255:         enddo 
256:  
257: ! read in position restraints 
258:        read(30,*)  
259:        read(30,*) SBMPRN 
260:        write(*,*) SBMPRN, ' position restraints' 
261:        do I=1,SBMPRN 
262:             read(30,*) SBMPRI(I),SBMPRK(I,1),SBMPRK(I,2),SBMPRK(I,3), 
263:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6), 
264:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3) 
265:             if(TARR(SBMPRI(I)) .ne. 0)then 
266:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I) 
267:                 call abort 
268:             endif 
269:             TARR(SBMPRI(I))=1 
270:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then   
271:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported' 
272:         endif 
273:        enddo  
274:  
275:        close(30)265:        close(30)
276:        end266:        end
277: 267: 
278: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^268: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
279: 269: 
280: 270: 
281: C271: C
282: C Calculate the Forces and energies272: C Calculate the Forces and energies
283: C273: C
284:       subroutine calc_energy_SBM(qo,natoms,GRAD, energy,Ib1, Ib2,274:       subroutine calc_energy_SBM(qo,natoms,GRAD, energy,Ib1, Ib2,
285:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,275:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,
286:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC, 276:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC, 
287:      Q NNCinc,NBA, NTA, NPA, NC,SBMCHARGE, SBMCHARGEON,NCswitch,277:      Q NNCsigma,NBA, NTA, NPA, NC,NCswitch,STT,NCcut,CONTACTTYPE,used)
288:      Q STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 
289:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 
290: 278: 
291:       INTEGER I, J, NATOMS,NBA, NTA, NPA, NC279:       INTEGER I, J, NATOMS,NBA, NTA, NPA, NC
292: 280: 
293:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY281:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
294:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)282:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
295: 283: 
296:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 284:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
297:      Q PHISBM(NPA), Sigma(NC),NCswitch,NCcut,STT, 285:      Q PHISBM(NPA), Sigma(NC),NCswitch,NCcut,STT, 
298:      Q EpsC(NC),SBMCHARGE(NATOMS)286:      Q EpsC(NC),  NNCsigma(NATOMS,NATOMS)
299:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 287:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 
300:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),SBMCHARGEON(NATOMS),288:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),
301:      Q PHITYPE(NPA),CONTACTTYPE,NNCinc(NATOMS,NATOMS)289:      Q PHITYPE(NPA),CONTACTTYPE
302:       DOUBLE PRECISION dx,dy,dz290:       DOUBLE PRECISION dx,dy,dz
303:       INTEGER SBMPRN, SBMPRI(NATOMS)291:         logical used(NATOMS,NATOMS)
304:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3) 
305:       do i = 1, natoms292:       do i = 1, natoms
306:         j = (i-1)*3293:          j = (i-1)*3
307:         x(i) = qo(j+1)294:          x(i) = qo(j+1)
308:         y(i) = qo(j+2)295:          y(i) = qo(j+2)
309:         z(i) = qo(j+3)296:          z(i) = qo(j+3)
310:         grad(j+1) = 0.0297:          grad(j+1) = 0.0
311:         grad(j+2) = 0.0298:         grad(j+2) = 0.0
312:         grad(j+3) = 0.0299:         grad(j+3) = 0.0
313:       enddo300:       enddo
314: 301: 
315:       energy = 0.0302:       energy = 0.0
316:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)303:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
317:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)304:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
318:         call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,305:         call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
319:      Q PHITYPE,PHISBM,NPA)306:      Q PHITYPE,PHISBM,NPA)
320:         call SBMContacts(x,y,z,grad, energy, natoms, IC, JC, 307:         call SBMContacts(x,y,z,grad, energy, natoms, IC, JC, 
321:      Q Sigma, EpsC, NC,CONTACTTYPE)308:      Q Sigma, EpsC, NC,CONTACTTYPE)
322:         call SBMNonContacts(x,y,z,grad, energy, natoms, 309:         call SBMNonContacts(x,y,z,grad, energy, natoms, 
323:      Q NNCinc,NCswitch,NCcut,STT)310:      Q NNCsigma,NCswitch,NCcut,STT,used)
324:         call SBMDHELEC(x,y,z,grad, energy, natoms,  
325:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON) 
326:         call SBMPR(x,y,z,grad, energy, natoms,  
327:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 
328: 311: 
329:       end312:       end
330: 313: 
331: 314: 
332: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<315: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
333: !* SBMBonds  computes the hookean force and energy between chosen atoms *316: !* SBMBonds  computes the hookean force and energy between chosen atoms *
334: !***********************************************************************317: !***********************************************************************
335: 318: 
336:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)319:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)
337:       USE KEY320:       USE KEY
338:       implicit NONE321:       implicit NONE
339:       integer I2, J2,  outE,I, N, J, NATOMS, NBA322:       integer I2, J2,  outE,I, N, J, NATOMS, NBA
340:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),323:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),
341:      Q energy324:      Q energy
342:       DOUBLE PRECISION r2, f, r1325:       DOUBLE PRECISION r2, f, r1
343:       DOUBLE PRECISION dx,dy,dz326:       DOUBLE PRECISION dx,dy,dz
344:       real f2327: 
345:         DOUBLE PRECISION Rb(NBA), bK(NBA)328:         DOUBLE PRECISION Rb(NBA), bK(NBA)
346:         INTEGER Ib1(NBA), Ib2(NBA)329:         INTEGER Ib1(NBA), Ib2(NBA)
347: 330: 
348: 331: 
349:         do 1 i=1, nBA332:         do 1 i=1, nBA
350:            I2 = Ib1(i)333:            I2 = Ib1(i)
351:            J2 = Ib2(i)334:            J2 = Ib2(i)
352:         dx = X(I2) - X(J2)335:         dx = X(I2) - X(J2)
353:         dy = Y(I2) - Y(J2)336:         dy = Y(I2) - Y(J2)
354:         dz = Z(I2) - Z(J2)337:         dz = Z(I2) - Z(J2)
355: 338: 
356:           r2 = dx**2 + dy**2 + dz**2339:           r2 = dx**2 + dy**2 + dz**2
357:           r1 = sqrt(r2)340:           r1 = sqrt(r2)
 341: 
358: ! energy calculation342: ! energy calculation
359:              Energy = Energy + bk(i)*(r1-Rb(i))**2/2.0343:              Energy = Energy + bk(i)*(r1-Rb(i))**2/2.0
360: 344: 
361: ! End energy calculation345: ! End energy calculation
362: 346: 
363: ! f_over_r is the force over the magnitude of r so there is no need to resolve347: ! f_over_r is the force over the magnitude of r so there is no need to resolve
364: ! the dx, dy and dz into unit vectors348: ! the dx, dy and dz into unit vectors
365: 349: 
366: ! the index i indicates the interaction between particle i and i+1350: ! the index i indicates the interaction between particle i and i+1
367: 351: 
368:              f = -bk(i)*(r1-Rb(i))/r1352:              f = -bk(i)*(r1-Rb(i))/r1
369:              !   f2=f*r1 
370: !            f = Rb(i)*bK(i)/r1 - bK(i)353: !            f = Rb(i)*bK(i)/r1 - bK(i)
371:         !write(*,*) i, f354:         !write(*,*) i, f
372:             ! now add the force355:             ! now add the force
373:               grad(I2*3-2) = grad(I2*3-2) - f * dx356:               grad(I2*3-2) = grad(I2*3-2) - f * dx
374:               grad(I2*3-1) = grad(I2*3-1) - f * dy357:               grad(I2*3-1) = grad(I2*3-1) - f * dy
375:               grad(I2*3)   = grad(I2*3)   - f * dz358:               grad(I2*3)   = grad(I2*3)   - f * dz
376: ! the negative sign is due to the computation of dx, dy and dz359: ! the negative sign is due to the computation of dx, dy and dz
377:               grad(J2*3-2) = grad(J2*3-2) + f * dx360:               grad(J2*3-2) = grad(J2*3-2) + f * dx
378:               grad(J2*3-1) = grad(J2*3-1) + f * dy361:               grad(J2*3-1) = grad(J2*3-1) + f * dy
379:               grad(J2*3)   = grad(J2*3)   + f * dz362:               grad(J2*3)   = grad(J2*3)   + f * dz
405:      + ZKJ, DF388:      + ZKJ, DF
406:       dimension  XIJ(NTA),YIJ(NTA),ZIJ(NTA),XKJ(NTA),YKJ(NTA),389:       dimension  XIJ(NTA),YIJ(NTA),ZIJ(NTA),XKJ(NTA),YKJ(NTA),
407:      + ZKJ(NTA),CST(NTA),EAW(NTA),RIJ(NTA),RKJ(NTA),RIK(NTA),390:      + ZKJ(NTA),CST(NTA),EAW(NTA),RIJ(NTA),RKJ(NTA),RIK(NTA),
408:      + DFW(NTA),ANT(NTA)391:      + DFW(NTA),ANT(NTA)
409:       DOUBLE PRECISION CT0, CT1, CT2, RIJ0, RKJ0, RIK0, ANT0, DA, ST, 392:       DOUBLE PRECISION CT0, CT1, CT2, RIJ0, RKJ0, RIK0, ANT0, DA, ST, 
410:      + CIK, CII, CKK, DT1, DT2, DT3, DT4, DT5, DT6, DT7, DT8, DT9, pt999393:      + CIK, CII, CKK, DT1, DT2, DT3, DT4, DT5, DT6, DT7, DT8, DT9, pt999
411:      Q , ebal,STH394:      Q , ebal,STH
412: 395: 
413: 396: 
414:         DOUBLE PRECISION ANTC(NTA), Tk(NTA)397:         DOUBLE PRECISION ANTC(NTA), Tk(NTA)
415:         INTEGER II, IT(NTA), JT(NTA), KT(NTA)398:         INTEGER JN, IT(NTA), JT(NTA), KT(NTA)
416:         INTEGER I3, J3, K3399:         INTEGER I3, J3, K3
417: 400: 
418:       data pt999 /1.0d0/401:       data pt999 /1.0d0/
419:       ebal= 0.0d0402:       ebal= 0.0d0
420: 403: 
421: 404: 
422:           DO II = 1, nTA405:           DO JN = 1, nTA
423:             I3 = IT(II)406:             I3 = IT(JN)
424:             J3 = JT(II)407:             J3 = JT(JN)
425:             K3 = KT(II)408:             K3 = KT(JN)
426: 409: 
427:             XIJ(II) = X(I3)-X(J3)410:             XIJ(JN) = X(I3)-X(J3)
428:             YIJ(II) = Y(I3)-Y(J3)411:             YIJ(JN) = Y(I3)-Y(J3)
429:             ZIJ(II) = Z(I3)-Z(J3)412:             ZIJ(JN) = Z(I3)-Z(J3)
430:             XKJ(II) = X(K3)-X(J3)413:             XKJ(JN) = X(K3)-X(J3)
431:             YKJ(II) = Y(K3)-Y(J3)414:             YKJ(JN) = Y(K3)-Y(J3)
432:             ZKJ(II) = Z(K3)-Z(J3)415:             ZKJ(JN) = Z(K3)-Z(J3)
433:           END DO416:           END DO
434:           DO II = 1,nTA417:           DO JN = 1,nTA
435:             RIJ0 = XIJ(II)*XIJ(II)+YIJ(II)*YIJ(II)+ZIJ(II)*ZIJ(II)418:             RIJ0 = XIJ(JN)*XIJ(JN)+YIJ(JN)*YIJ(JN)+ZIJ(JN)*ZIJ(JN)
436:             RKJ0 = XKJ(II)*XKJ(II)+YKJ(II)*YKJ(II)+ZKJ(II)*ZKJ(II)419:             RKJ0 = XKJ(JN)*XKJ(JN)+YKJ(JN)*YKJ(JN)+ZKJ(JN)*ZKJ(JN)
437:             RIK0 = SQRT(RIJ0*RKJ0)420:             RIK0 = SQRT(RIJ0*RKJ0)
438:             CT0 = (XIJ(II)*XKJ(II)+YIJ(II)*YKJ(II)+ZIJ(II)*ZKJ(II))/RIK0421:             CT0 = (XIJ(JN)*XKJ(JN)+YIJ(JN)*YKJ(JN)+ZIJ(JN)*ZKJ(JN))/RIK0
439:             CT1 = MAX(-pt999,CT0)422:             CT1 = MAX(-pt999,CT0)
440:             CT2 = MIN(pt999,CT1)423:             CT2 = MIN(pt999,CT1)
441:             CST(II) = CT2424:             CST(JN) = CT2
442:             ANT(II) = ACOS(CT2)425:             ANT(JN) = ACOS(CT2)
443:             RIJ(II) = RIJ0426:             RIJ(JN) = RIJ0
444:             RKJ(II) = RKJ0427:             RKJ(JN) = RKJ0
445:             RIK(II) = RIK0428:             RIK(JN) = RIK0
446:           END DO429:           END DO
447: 430: 
448: ! end of insertion431: ! end of insertion
449: 432: 
450: 433: 
451:           DO II = 1,nTA434:           DO JN = 1,nTA
452:             ANT0 = ANT(II)435:             ANT0 = ANT(JN)
453:             DA = ANT0 - ANTC(II)436:             DA = ANT0 - ANTC(JN)
454:             DF = TK(II)*DA437:             DF = TK(JN)*DA
455: 438: 
456: 439: 
457:             DFW(II) = -(DF)/SIN(ANT0)440:             DFW(JN) = -(DF)/SIN(ANT0)
458: 441: 
459:           END DO442:           END DO
460:           DO II = 1,nTA443:           DO JN = 1,nTA
461:             I3 = IT(II)444:             I3 = IT(JN)
462:             J3 = JT(II)445:             J3 = JT(JN)
463:             K3 = KT(II)446:             K3 = KT(JN)
464:             ST = DFW(II)447:             ST = DFW(JN)
465:             STH = ST*CST(II)448:             STH = ST*CST(JN)
466:             CIK = ST/RIK(II)449:             CIK = ST/RIK(JN)
467:             CII = STH/RIJ(II)450:             CII = STH/RIJ(JN)
468:             CKK = STH/RKJ(II)451:             CKK = STH/RKJ(JN)
469:             DT1 = CIK*XKJ(II)-CII*XIJ(II)452:             DT1 = CIK*XKJ(JN)-CII*XIJ(JN)
470:             DT2 = CIK*YKJ(II)-CII*YIJ(II)453:             DT2 = CIK*YKJ(JN)-CII*YIJ(JN)
471:             DT3 = CIK*ZKJ(II)-CII*ZIJ(II)454:             DT3 = CIK*ZKJ(JN)-CII*ZIJ(JN)
472:             DT7 = CIK*XIJ(II)-CKK*XKJ(II)455:             DT7 = CIK*XIJ(JN)-CKK*XKJ(JN)
473:             DT8 = CIK*YIJ(II)-CKK*YKJ(II)456:             DT8 = CIK*YIJ(JN)-CKK*YKJ(JN)
474:             DT9 = CIK*ZIJ(II)-CKK*ZKJ(II)457:             DT9 = CIK*ZIJ(JN)-CKK*ZKJ(JN)
475:             DT4 = -DT1-DT7458:             DT4 = -DT1-DT7
476:             DT5 = -DT2-DT8459:             DT5 = -DT2-DT8
477:             DT6 = -DT3-DT9460:             DT6 = -DT3-DT9
478: C461: C
479: 462: 
480:             grad(I3*3-2) = grad(I3*3-2)+ DT1463:             grad(I3*3-2) = grad(I3*3-2)+ DT1
481:             grad(I3*3-1) = grad(I3*3-1)+ DT2464:             grad(I3*3-1) = grad(I3*3-1)+ DT2
482:             grad(I3*3)   = grad(I3*3)  + DT3465:             grad(I3*3)   = grad(I3*3)  + DT3
483:             grad(J3*3-2) = grad(J3*3-2)+ DT4466:             grad(J3*3-2) = grad(J3*3-2)+ DT4
484:             grad(J3*3-1) = grad(J3*3-1)+ DT5467:             grad(J3*3-1) = grad(J3*3-1)+ DT5
503: 486: 
504: 487: 
505: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<488: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
506: !* SBMdihedral computes the dihedral angles and the forces due to them *489: !* SBMdihedral computes the dihedral angles and the forces due to them *
507: !**********************************************************************490: !**********************************************************************
508: 491: 
509:       SUBROUTINE SBMdihedral(x,y,z,grad, energy, NATOMS,IP,JP,KP,LP,PK,492:       SUBROUTINE SBMdihedral(x,y,z,grad, energy, NATOMS,IP,JP,KP,LP,PK,
510:      Q PHITYPE,PHISBM,NPA)493:      Q PHITYPE,PHISBM,NPA)
511:       USE KEY494:       USE KEY
512:       implicit NONE495:       implicit NONE
513:       integer I, N, J, NATOMS, NPA, II496:       integer I, N, J, NATOMS, NPA, JN
514:       DOUBLE PRECISION x(NATOMS),y(NATOMS),z(NATOMS),497:       DOUBLE PRECISION x(NATOMS),y(NATOMS),z(NATOMS),
515:      Q grad(3*NATOMS),energy498:      Q grad(3*NATOMS),energy
516: 499: 
517: 500: 
518:       DOUBLE PRECISION PK(NPA),PHISBM(NPA)501:       DOUBLE PRECISION PK(NPA),PHISBM(NPA)
519:       INTEGER IP(NPA), JP(NPA), KP(NPA), LP(NPA) 502:       INTEGER IP(NPA), JP(NPA), KP(NPA), LP(NPA) 
520:       INTEGER PHITYPE(NPA)503:       INTEGER PHITYPE(NPA)
521: 504: 
522:       double precision lfac505:       double precision lfac
523:       integer I3, J3, K3, L3506:       integer I3, J3, K3, L3
524:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,507:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,
525:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,508:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,
526:      + DZ, GX,GY,GZ,CT,CPHI,509:      + DZ, GX,GY,GZ,CT,CPHI,
527:      + SPHI,Z1, Z2,FXI,FYI,FZI,510:      + SPHI,Z1, Z2,FXI,FYI,FZI,
528:      + FXJ,FYJ,FZJ, FXK,FYK,FZK,511:      + FXJ,FYJ,FZJ, FXK,FYK,FZK,
529:      + FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,ftem,CT0,CT1,AP0,AP1,512:      + FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,ftem,CT0,CT1,AP0,AP1,
530:      + Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,513:      + Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,
531:      + S,HGoverG,FGoverG,A1,A3514:      +  DC1, DC2, DC3, DC4, DC5, DC6,S,HGoverG,FGoverG,A1,A3
532: 515: 
533:       DIMENSION XIJ(NPA),YIJ(NPA),ZIJ(NPA),RIJ(NPA),XKJ(NPA),516:       DIMENSION XIJ(NPA),YIJ(NPA),ZIJ(NPA),RIJ(NPA),XKJ(NPA),
534:      + YKJ(NPA),517:      + YKJ(NPA),
535:      + ZKJ(NPA),RKJ(NPA),XKL(NPA),YKL(NPA),ZKL(NPA),RKL(NPA),DX(NPA),518:      + ZKJ(NPA),RKJ(NPA),XKL(NPA),YKL(NPA),ZKL(NPA),RKL(NPA),DX(NPA),
536:      + DY(NPA),519:      + DY(NPA),
537:      + DZ(NPA), GX(NPA),GY(NPA),GZ(NPA),CT(NPA),CPHI(NPA),520:      + DZ(NPA), GX(NPA),GY(NPA),GZ(NPA),CT(NPA),CPHI(NPA),
538:      + SPHI(NPA),Z1(NPA), Z2(NPA),FXI(NPA),FYI(NPA),FZI(NPA),521:      + SPHI(NPA),Z1(NPA), Z2(NPA),FXI(NPA),FYI(NPA),FZI(NPA),
539:      + FXJ(NPA),FYJ(NPA),FZJ(NPA), FXK(NPA),FYK(NPA),FZK(NPA),522:      + FXJ(NPA),FYJ(NPA),FZJ(NPA), FXK(NPA),FYK(NPA),FZK(NPA),
540:      + FXL(NPA),FYL(NPA),FZL(NPA),DF(NPA),FGoverG(NPA),HGoverG(NPA)523:      + FXL(NPA),FYL(NPA),FZL(NPA),DF(NPA),FGoverG(NPA),HGoverG(NPA)
541: C524: C
543: 526: 
544:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,527:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,
545:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z528:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z
546: 529: 
547:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/530:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/
548:       data zero,one,two,four,six,twelve/0.d0,1.d0,2.d0,4.d0,6.d0,12.d0/531:       data zero,one,two,four,six,twelve/0.d0,1.d0,2.d0,4.d0,6.d0,12.d0/
549: 532: 
550:       double precision pi533:       double precision pi
551:       pi = 3.14159265358979323846264338327950288419716939937510534:       pi = 3.14159265358979323846264338327950288419716939937510
552: 535: 
553:           DO II = 1,nPA536:           DO JN = 1,nPA
554: 537: 
555:             I3 = IP(II)538:             I3 = IP(JN)
556:             J3 = JP(II)539:             J3 = JP(JN)
557:             K3 = KP(II)540:             K3 = KP(JN)
558:             L3 = LP(II)541:             L3 = LP(JN)
559: 542: 
560:  543:  
561: 544: 
562:             XIJ(II) = X(I3)-X(J3)545:             XIJ(JN) = X(I3)-X(J3)
563:             YIJ(II) = Y(I3)-Y(J3)546:             YIJ(JN) = Y(I3)-Y(J3)
564:             ZIJ(II) = Z(I3)-Z(J3)547:             ZIJ(JN) = Z(I3)-Z(J3)
565:             XKJ(II) = X(K3)-X(J3)548:             XKJ(JN) = X(K3)-X(J3)
566:             YKJ(II) = Y(K3)-Y(J3)549:             YKJ(JN) = Y(K3)-Y(J3)
567:             ZKJ(II) = Z(K3)-Z(J3)550:             ZKJ(JN) = Z(K3)-Z(J3)
568:             RKJ(II) = sqrt(XKJ(II)**2+YKJ(II)**2+ZKJ(II)**2)551:             RKJ(JN) = sqrt(XKJ(JN)**2+YKJ(JN)**2+ZKJ(JN)**2)
569:             XKL(II) = X(K3)-X(L3)552:             XKL(JN) = X(K3)-X(L3)
570:             YKL(II) = Y(K3)-Y(L3)553:             YKL(JN) = Y(K3)-Y(L3)
571:             ZKL(II) = Z(K3)-Z(L3)                                  554:             ZKL(JN) = Z(K3)-Z(L3)                                  
572: 555: 
573: 556: 
574:             FGoverG(II)=-(XIJ(II)*XKJ(II)+YIJ(II)*YKJ(II)+557:             FGoverG(JN)=-(XIJ(JN)*XKJ(JN)+YIJ(JN)*YKJ(JN)+
575:      Q   ZIJ(II)*ZKJ(II))/RKJ(II)558:      Q   ZIJ(JN)*ZKJ(JN))/RKJ(JN)
576:             HGoverG(II)=(XKL(II)*XKJ(II)+YKL(II)*YKJ(II)+559:             HGoverG(JN)=(XKL(JN)*XKJ(JN)+YKL(JN)*YKJ(JN)+
577:      Q   ZKL(II)*ZKJ(II))/RKJ(II)560:      Q   ZKL(JN)*ZKJ(JN))/RKJ(JN)
578:           END DO561:           END DO
579: C DX is the M vector and G is the N vector562: C DX is the M vector and G is the N vector
580:           DO II = 1,nPA563:           DO JN = 1,nPA
581:             DX(II) = YIJ(II)*ZKJ(II)-ZIJ(II)*YKJ(II)564:             DX(JN) = YIJ(JN)*ZKJ(JN)-ZIJ(JN)*YKJ(JN)
582:             DY(II) = ZIJ(II)*XKJ(II)-XIJ(II)*ZKJ(II)565:             DY(JN) = ZIJ(JN)*XKJ(JN)-XIJ(JN)*ZKJ(JN)
583:             DZ(II) = XIJ(II)*YKJ(II)-YIJ(II)*XKJ(II)566:             DZ(JN) = XIJ(JN)*YKJ(JN)-YIJ(JN)*XKJ(JN)
584:             GX(II) = ZKJ(II)*YKL(II)-YKJ(II)*ZKL(II)567:             GX(JN) = ZKJ(JN)*YKL(JN)-YKJ(JN)*ZKL(JN)
585:             GY(II) = XKJ(II)*ZKL(II)-ZKJ(II)*XKL(II)568:             GY(JN) = XKJ(JN)*ZKL(JN)-ZKJ(JN)*XKL(JN)
586:             GZ(II) = YKJ(II)*XKL(II)-XKJ(II)*YKL(II)569:             GZ(JN) = YKJ(JN)*XKL(JN)-XKJ(JN)*YKL(JN)
587:           END DO570:           END DO
588: C571: C
589: 572: 
590: 573: 
591: ! so far so good574: ! so far so good
592: 575: 
593: 576: 
594:           DO II = 1,nPA577:           DO JN = 1,nPA
595:             FXI(II) = SQRT(DX(II)*DX(II)578:             FXI(JN) = SQRT(DX(JN)*DX(JN)
596:      Q                    +DY(II)*DY(II)579:      Q                    +DY(JN)*DY(JN)
597:      Q                    +DZ(II)*DZ(II))580:      Q                    +DZ(JN)*DZ(JN))
598:             FYI(II) = SQRT(GX(II)*GX(II)581:             FYI(JN) = SQRT(GX(JN)*GX(JN)
599:      Q                    +GY(II)*GY(II)582:      Q                    +GY(JN)*GY(JN)
600:      Q                    +GZ(II)*GZ(II))583:      Q                    +GZ(JN)*GZ(JN))
601:             CT(II) = DX(II)*GX(II)+DY(II)*GY(II)+DZ(II)*GZ(II)584:             CT(JN) = DX(JN)*GX(JN)+DY(JN)*GY(JN)+DZ(JN)*GZ(JN)
602:           END DO585:           END DO
603: 586: 
604:          DO II = 1,nPA587:          DO JN = 1,nPA
605:             z10 = 1.0/FXI(II)588:             z10 = 1.0/FXI(jn)
606:             z20 = 1.0/FYI(II)589:             z20 = 1.0/FYI(jn)
607:             Z12 = Z10*Z20590:             Z12 = Z10*Z20
608:             Z1(II) = Z10591:             Z1(JN) = Z10
609:             Z2(II) = Z20592:             Z2(JN) = Z20
610:             ftem = zero593:             ftem = zero
611:             CT0 = MIN(one,CT(II)*Z12)594:             CT0 = MIN(one,CT(JN)*Z12)
612:             CT1 = MAX(-one,CT0)595:             CT1 = MAX(-one,CT0)
613:             S = XKJ(II)*(DZ(II)*GY(II)-DY(II)*GZ(II))+596:             S = XKJ(JN)*(DZ(JN)*GY(JN)-DY(JN)*GZ(JN))+
614:      Q          YKJ(II)*(DX(II)*GZ(II)-DZ(II)*GX(II))+597:      Q          YKJ(JN)*(DX(JN)*GZ(JN)-DZ(JN)*GX(JN))+
615:      Q          ZKJ(II)*(DY(II)*GX(II)-DX(II)*GY(II))598:      Q          ZKJ(JN)*(DY(JN)*GX(JN)-DX(JN)*GY(JN))
616:             AP0 = ACOS(CT1)599:             AP0 = ACOS(CT1)
617:             AP1 = PI-SIGN(AP0,S)600:             AP1 = PI-SIGN(AP0,S)
618: 601: 
619:             CT(II) = AP1602:             CT(JN) = AP1
620:             CPHI(II) = COS(AP1)603:             CPHI(JN) = COS(AP1)
621:             SPHI(II) = SIN(AP1)604:             SPHI(JN) = SIN(AP1)
622:          END DO605:          END DO
623: 606: 
624: 607: 
625:         DO II = 1,nPA608:         DO JN = 1,nPA
626: !            CT0 = CT(II)609: !            CT0 = CT(JN)
627: ! Here is the energy part610: ! Here is the energy part
628:           A1=CT(II)-PHISBM(II)611:           A1=CT(JN)-PHISBM(JN)
629:           A3=A1*3612:           A3=A1*3
630: 613: 
631:         if(PHITYPE(II) .eq. 1)then614:         if(PHITYPE(JN) .eq. 1)then
632:           Energy =  Energy + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))615:           Energy =  Energy + PK(JN)*(3.0/2.0-cos(A1)-0.5*cos(A3))
633: ! dE/dPHI616: ! dE/dPHI
634:             DF(II)=PK(II)*(sin(A1)+1.5*sin(A3))617:             DF(JN)=PK(JN)*(sin(A1)+1.5*sin(A3))
635:         elseif(PHITYPE(II) .eq. 2)then618:         elseif(PHITYPE(JN) .eq. 2)then
636:           if(A1 .gt. PI)then619:           if(A1 .gt. PI)then
637:                 A1=A1-2*PI620:                 A1=A1-2*PI
638:           elseif(A1 .lt. -PI)then621:           elseif(A1 .lt. -PI)then
639:                 A1=A1+2*PI622:                 A1=A1+2*PI
640:           endif623:           endif
641: 624: 
642:           Energy =  Energy + 0.5*PK(II)*A1**2625:           Energy =  Energy + PK(JN)*A1**2
643: ! dE/dPHI626: ! dE/dPHI
644:             DF(II)=PK(II)*A1627:             DF(JN)=2*PK(JN)*A1
645: 628: 
646:         else629:         else
647: 630: 
648:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)631:         write(*,*) 'unrecognized dihedral type', PHITYPE(JN)
649:         STOP632:         STOP
650:         endif633:         endif
651:         END DO634:         END DO
652: 635: 
653: 636: 
654: 637: 
655: ! insert the new 638: ! insert the new 
656: 639: 
657:        ! now, do dPhi/dX640:        ! now, do dPhi/dX
658: 641: 
659:          DO II = 1,nPA642:          DO JN = 1,nPA
660: 643: 
661: ! |G|/|A|**2 644: ! |G|/|A|**2 
662:             TT1 = Z1(II)*Z1(II)*RKJ(II)*DF(II)645:             TT1 = Z1(JN)*Z1(JN)*RKJ(JN)*DF(JN)
663: ! FG/(A**2*|G|)646: ! FG/(A**2*|G|)
664:             TT2 = FGoverG(II)*Z1(II)*Z1(II)*DF(II)647:             TT2 = FGoverG(JN)*Z1(JN)*Z1(JN)*DF(JN)
665: ! HG/(B**2*|G|)648: ! HG/(B**2*|G|)
666:             TT3 = HGoverG(II)*Z2(II)*Z2(II)*DF(II)649:             TT3 = HGoverG(JN)*Z2(JN)*Z2(JN)*DF(JN)
667: ! |G|/|B|**2 650: ! |G|/|B|**2 
668:             TT4 = Z2(II)*Z2(II)*RKJ(II)*DF(II)651:             TT4 = Z2(JN)*Z2(JN)*RKJ(JN)*DF(JN)
669: 652: 
670: 653: 
671: ! note: negatives are flipped from paper because A=-DX654: ! note: negatives are flipped from paper because A=-DX
672:         TT1X=TT1*DX(II)655:         TT1X=TT1*DX(JN)
673:         TT1Y=TT1*DY(II)656:         TT1Y=TT1*DY(JN)
674:         TT1Z=TT1*DZ(II)657:         TT1Z=TT1*DZ(JN)
675: 658: 
676:         TT2X=TT2*DX(II)659:         TT2X=TT2*DX(JN)
677:         TT2Y=TT2*DY(II)660:         TT2Y=TT2*DY(JN)
678:         TT2Z=TT2*DZ(II)661:         TT2Z=TT2*DZ(JN)
679: 662: 
680: 663: 
681:         TT3X=TT3*GX(II)664:         TT3X=TT3*GX(JN)
682:         TT3Y=TT3*GY(II)665:         TT3Y=TT3*GY(JN)
683:         TT3Z=TT3*GZ(II)666:         TT3Z=TT3*GZ(JN)
684: 667: 
685: 668: 
686:         TT4X=TT4*GX(II)669:         TT4X=TT4*GX(JN)
687:         TT4Y=TT4*GY(II)670:         TT4Y=TT4*GY(JN)
688:         TT4Z=TT4*GZ(II)671:         TT4Z=TT4*GZ(JN)
689: 672: 
690:             I3 = IP(II)673:             I3 = IP(JN)
691:             J3 = JP(II)674:             J3 = JP(JN)
692:             K3 = KP(II)675:             K3 = KP(JN)
693:             L3 = LP(II)676:             L3 = LP(JN)
694: 677: 
695: 678: 
696:             grad(I3*3-2) =  grad(I3*3-2)  + TT1X  679:             grad(I3*3-2) =  grad(I3*3-2)  + TT1X  
697:             grad(I3*3-1) =  grad(I3*3-1)  + TT1Y680:             grad(I3*3-1) =  grad(I3*3-1)  + TT1Y
698:             grad(I3*3)   =  grad(I3*3)    + TT1Z681:             grad(I3*3)   =  grad(I3*3)    + TT1Z
699:             grad(J3*3-2) =  grad(J3*3-2)  - TT1X - TT2X - TT3X682:             grad(J3*3-2) =  grad(J3*3-2)  - TT1X - TT2X - TT3X
700:             grad(J3*3-1) =  grad(J3*3-1)  - TT1Y - TT2Y - TT3Y683:             grad(J3*3-1) =  grad(J3*3-1)  - TT1Y - TT2Y - TT3Y
701:             grad(J3*3)   =  grad(J3*3)    - TT1Z - TT2Z - TT3Z684:             grad(J3*3)   =  grad(J3*3)    - TT1Z - TT2Z - TT3Z
702:             grad(K3*3-2) =  grad(K3*3-2)  + TT2X + TT3X - TT4X685:             grad(K3*3-2) =  grad(K3*3-2)  + TT2X + TT3X - TT4X
703:             grad(K3*3-1) =  grad(K3*3-1)  + TT2Y + TT3Y - TT4Y686:             grad(K3*3-1) =  grad(K3*3-1)  + TT2Y + TT3Y - TT4Y
729:      Q , grad(3*NATOMS), energy712:      Q , grad(3*NATOMS), energy
730:       DOUBLE PRECISION dx,dy,dz713:       DOUBLE PRECISION dx,dy,dz
731: 714: 
732:       integer C1, C2, ConfID, Q, SC1, SC2, Cf1, cf2715:       integer C1, C2, ConfID, Q, SC1, SC2, Cf1, cf2
733:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r, dsig, deps, 716:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r, dsig, deps, 
734:      Q s1, s2, ep1, ep2, r1, rc,r, summm717:      Q s1, s2, ep1, ep2, r1, rc,r, summm
735: 718: 
736:         DOUBLE PRECISION Sigma(NC), EpsC(NC)719:         DOUBLE PRECISION Sigma(NC), EpsC(NC)
737:         INTEGER IC(NC), JC(NC),CONTACTTYPE720:         INTEGER IC(NC), JC(NC),CONTACTTYPE
738: 721: 
739: ! implement type 3 as gaussian 
740: ! implement type 4 as dual-basin gaussian 
741: ! question: Should type be part of the interaction, or should they be organized 
742: ! into separate lists? 
743:  
744: 722: 
745: ! type 1 is 6-12 interaction723: ! type 1 is 6-12 interaction
746:         if(CONTACTTYPE .eq. 1)then724:         if(CONTACTTYPE .eq. 1)then
747: 725: 
748:                 do i=1, NC726:        do i=1, NC
749:                727:        
750:                         C1 = IC(i)728:         C1 = IC(i)
751:                         C2 = JC(i)729:         C2 = JC(i)
752:                         dx = X(C1) - X(C2)730:         dx = X(C1) - X(C2)
753:                         dy = Y(C1) - Y(C2)731:         dy = Y(C1) - Y(C2)
754:                         dz = Z(C1) - Z(C2)732:         dz = Z(C1) - Z(C2)
755:                 733: 
756:                         r2 = dx**2 + dy**2 + dz**2734:           r2 = dx**2 + dy**2 + dz**2
757:                         rm2 = 1.0/r2735: 
758:                         rm2 = rm2*sigma(i)736:               rm2 = 1.0/r2
759:                         rm6 = rm2**3737:               rm2 = rm2*sigma(i)
760:                 738:               rm6 = rm2**3
761:                         energy = energy + epsC(i)*rm6*(rm6-2.0)739: 
762:                 740:         energy = energy + epsC(i)*rm6*(rm6-2.0)
763:                         f_over_r = -epsC(i)*12.0*rm6*(rm6-1.0)/r2741: 
764:                 742:         f_over_r = -epsC(i)*12.0*rm6*(rm6-1.0)/r2
765:                         grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx743: 
766:                         grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy744:               grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
767:                         grad(3*C1)   = grad(3*C1)   + f_over_r * dz745:               grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
768:                 746:               grad(3*C1)   = grad(3*C1)   + f_over_r * dz
769:                         grad(3*C2-2) =  grad(3*C2-2) - f_over_r * dx747: 
770:                         grad(3*C2-1) =  grad(3*C2-1) - f_over_r * dy748:               grad(3*C2-2) =  grad(3*C2-2) - f_over_r * dx
771:                         grad(3*C2)   =  grad(3*C2)   - f_over_r * dz749:               grad(3*C2-1) =  grad(3*C2-1) - f_over_r * dy
772:                 enddo750:               grad(3*C2)   =  grad(3*C2)   - f_over_r * dz
 751:               enddo
773: 752: 
774: ! type 2 is 10-12 interaction753: ! type 2 is 10-12 interaction
775:         elseif(CONTACTTYPE .eq. 2)then754:         elseif(CONTACTTYPE .eq. 2)then
776:         do i=1, NC755:        do i=1, NC
777:        756:        
778:                 C1 = IC(i)757:         C1 = IC(i)
779:                 C2 = JC(i)758:         C2 = JC(i)
780:                 dx = X(C1) - X(C2)759:         dx = X(C1) - X(C2)
781:                 dy = Y(C1) - Y(C2)760:         dy = Y(C1) - Y(C2)
782:                 dz = Z(C1) - Z(C2)761:         dz = Z(C1) - Z(C2)
783:         762: 
784:                 r2 = dx**2 + dy**2 + dz**2763:           r2 = dx**2 + dy**2 + dz**2
785:         764: 
786:                 rm2 = 1.0/r2765:               rm2 = 1.0/r2
787:                 rm2 = rm2*sigma(i)766:               rm2 = rm2*sigma(i)
788:                 rm10 = rm2**5767:               rm10 = rm2**5
789:                 768: 
790:                 energy = energy + epsC(i)*rm10*(5*rm2-6.0)769:         energy = energy + epsC(i)*rm10*(5*rm2-6.0)
791:                 f_over_r = -epsc(i)*60.0*rm10*(rm2-1.0)/r2770:         f_over_r = -epsc(i)*60.0*rm10*(rm2-1.0)/r2
792:                 771: 
793:                 grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx772:               grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
794:                 grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy773:               grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
795:                 grad(3*C1)   = grad(3*C1)   + f_over_r * dz774:               grad(3*C1)   = grad(3*C1)   + f_over_r * dz
796:                 775: 
797:                 grad(3*C2-2) =  grad(3*C2-2) - f_over_r * dx776:               grad(3*C2-2) =  grad(3*C2-2) - f_over_r * dx
798:                 grad(3*C2-1) =  grad(3*C2-1) - f_over_r * dy777:               grad(3*C2-1) =  grad(3*C2-1) - f_over_r * dy
799:                 grad(3*C2)   =  grad(3*C2)   - f_over_r * dz778:               grad(3*C2)   =  grad(3*C2)   - f_over_r * dz
800:               enddo779:               enddo
801: 780: 
802:         else781:         else
803:         write(*,*) CONTACTTYPE, 'is not a valid contact selection'782:         write(*,*) CONTACTTYPE, 'is not a valid contact selection'
804:         stop783:         stop
805:         endif784:         endif
806: 785: 
807: 786: 
808:       end787:       end
809: 788: 
810: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^789: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
811: 790: 
812: 791: 
813: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<792: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
814: !* SBMNonContacts computes the forces due to non native contacts       *793: !* SBMNonContacts computes the forces due to non native contacts       *
815: !**********************************************************************794: !**********************************************************************
816: 795: 
817:       subroutine SBMnoncontacts(x,y,z,grad, energy, 796:       subroutine SBMnoncontacts(x,y,z,grad, energy, 
818:      Q NATOMS,NNCinc,NCswitch,NCcut,STT)797:      Q NATOMS,NNCsigma,NCswitch,NCcut,STT,used)
819:       USE KEY798:       USE KEY
820:       implicit NONE799:       implicit NONE
821:       integer I, N, J, AN, NATOMS800:       integer I, N, J, AN, NATOMS
822: 801: 
823: 802: 
824:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),803:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
825:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC804:      Q grad(3*NATOMS), energy,STT,SA,SB,SC
826: 805: 
827:       integer C1, C2, ii,jj,kk,k,l,iii,jjj806:       integer C1, C2, ii,jj,kk,k,l,iii,jjj
828:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 807:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
829: 808: 
830:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS, 
831:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS 
832: 809: 
833:         INTEGER NNCinc(NATOMS,NATOMS)810:         DOUBLE PRECISION NNCsigma(NATOMS,NATOMS)
834:         DOUBLE PRECISION dx,dy,dz811:       DOUBLE PRECISION dx,dy,dz
835:         integer tempN, alpha812:         integer tempN, temparray,alpha
 813:         logical used
 814:         dimension used(NATOMS,NATOMS)
836:         double precision Rdiff,Vfunc,Ffunc815:         double precision Rdiff,Vfunc,Ffunc
837:         double precision Rcut2,Rswitch2816:         double precision Rcut2,Rswitch2
838:         ! Ngrid is the number of atoms in that grid point817:         ! Ngrid is the number of atoms in that grid point
839:         ! grid is the array of atoms in each grid818:         ! grid is the array of atoms in each grid
840:         integer Ngrid,grid,maxgrid,maxpergrid819:         integer Ngrid,grid,maxgrid,maxpergrid
841:         ! number of atoms per grid, max820:         ! number of atoms per grid, max
842:         parameter (maxpergrid=50)821:         parameter (maxpergrid=2000)
 822:         dimension temparray(NATOMS*NATOMS,2)
843:         ! dimensions of grid823:         ! dimensions of grid
844:         parameter (maxgrid=15)824:         parameter (maxgrid=15)
845:         dimension Ngrid(maxgrid,maxgrid,maxgrid),825:         dimension Ngrid(maxgrid,maxgrid,maxpergrid),
846:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)826:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
847:         integer MaxGridX,MaxGridY,MaxGridZ827:         integer MaxGridX,MaxGridY,MaxGridZ
848:         double precision gridsize,RD1828:         double precision gridsize,RD1
849:         double precision minX,minY,minZ,maxX,maxY,maxZ829:         double precision minX,minY,minZ,maxX,maxY,maxZ
850:         integer Xgrid,Ygrid,Zgrid,TID830:         integer Xgrid,Ygrid,Zgrid
851:         double precision Etemp831:         double precision Etemp
852:         Rdiff=NCcut-NCswitch832:         Rdiff=NCcut-NCswitch
853:         alpha=12833:         alpha=12
854: 834: 
855: 835: 
856:         GRIDSIZE=NCcut*1.01836:         GRIDSIZE=NCcut*1.01
857:         Rcut2=NCcut**2837:         Rcut2=NCcut**2
858:         Rswitch2=NCswitch**2838:         Rswitch2=NCswitch**2
859: 839: 
860:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 840:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
861:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )841:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
862: 842: 
863:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)843:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)
864: 844: 
865:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)845:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)
866: 846: 
867: 847: 
868: 848: 
869: ! grid the system849: 
 850: !! make a neighbor list
 851: 
870:         minX=10000000852:         minX=10000000
871:         minY=10000000853:         minY=10000000
872:         minZ=10000000854:         minZ=10000000
873:         855:         
874:         maxX=-10000000856:         maxX=-10000000
875:         maxY=-10000000857:         maxY=-10000000
876:         maxZ=-10000000858:         maxZ=-10000000
877: 859: 
878:         do i=1,NATOMS860:         do i=1,NATOMS
879:         gradOMP(i*3)=0 
880:         gradOMP(i*3-1)=0 
881:         gradOMP(i*3-2)=0 
882:  
883:            if(X(i) .gt. maxX)then861:            if(X(i) .gt. maxX)then
884:                 maxX=X(i)862:                 maxX=X(i)
885:            elseif(Y(i) .gt. maxY)then863:            endif
 864: 
 865:            if(Y(i) .gt. maxY)then
886:                 maxY=Y(i)866:                 maxY=Y(i)
887:            endif867:            endif
888: 868: 
889:            if(Z(i) .gt. maxZ)then869:            if(Z(i) .gt. maxZ)then
890:                 maxZ=Z(i)870:                 maxZ=Z(i)
891:            elseif(X(i) .lt. minX)then871:            endif
 872: 
 873:            if(X(i) .lt. minX)then
892:                 minX=X(i)874:                 minX=X(i)
893:            endif875:            endif
894: 876: 
895:            if(Y(i) .lt. minY)then877:            if(Y(i) .lt. minY)then
896:                 minY=Y(i)878:                 minY=Y(i)
897:            elseif(Z(i) .lt. minZ)then879:            endif
 880: 
 881:            if(Z(i) .lt. minZ)then
898:                 minZ=Z(i)882:                 minZ=Z(i)
899:            endif883:            endif
900:         enddo884:         enddo
901: 885: 
902:         maxgridX=int((maxX-minX)/gridsize)+1886:         maxgridX=int((maxX-minX)/gridsize)+1
903:         maxgridY=int((maxY-minY)/gridsize)+1887:         maxgridY=int((maxY-minY)/gridsize)+1
904:         maxgridZ=int((maxZ-minZ)/gridsize)+1888:         maxgridZ=int((maxZ-minZ)/gridsize)+1
 889: !        write(*,*) gridsize, maxgrid
 890: !        write(*,*) maxgridX,maxgridY,maxgridZ
 891: !        write(*,*) minX,maxX
 892: !        write(*,*) minY,maxY
 893: !        write(*,*) minZ,maxZ
905: 894: 
906:         if(maxgridX .ge. maxgrid .or. 895:         if(maxgridX .ge. maxgrid .or. 
907:      Q  maxgridY .ge. maxgrid .or.896:      Q  maxgridY .ge. maxgrid .or.
908:      Q  maxgridZ .ge. maxgrid )then897:      Q  maxgridZ .ge. maxgrid )then
909:         write(*,*) 'system got too big for grid searching...'898:         write(*,*) 'system got too big for grid searching...'
910:         call abort899: !call abort
911:         endif900:         endif
912: 901: 
913:  
914: !!  Add a second grid that only includes atoms that are charged 
915:  
916:         do i=1,maxgrid902:         do i=1,maxgrid
917:          do j=1,maxgrid903:          do j=1,maxgrid
918:           do k=1,maxgrid904:           do k=1,maxgrid
919:                 Ngrid(i,j,k)=0905:                 Ngrid(i,j,k)=0
920:           enddo906:           enddo
921:          enddo907:          enddo
922:         enddo908:         enddo
923:         do i=1,NATOMS909:         do i=1,NATOMS
924: 910: 
925:                 Xgrid=int((X(i)-minX)/gridsize)+1911:                 Xgrid=int((X(i)-minX)/gridsize)+1
926:                 Ygrid=int((Y(i)-minY)/gridsize)+1912:                 Ygrid=int((Y(i)-minY)/gridsize)+1
927:                 Zgrid=int((Z(i)-minZ)/gridsize)+1913:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
 914:         !        write(*,*) Xgrid,Ygrid,Zgrid
928:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1915:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
 916: !                write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
929:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then917:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
930:                         write(*,*) 'too many atoms in a grid'918:                         write(*,*) 'too many atoms in a grid'
931:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid919:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
932:                         call abort920: !                call abort
933:                 endif921:                 endif
934:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i922:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
935:         enddo923:         enddo
936: 924: 
937: 925: 
938:            tempN=0 
939:  
940: ! add a second loop that goes over charged atoms 
941:  
942: !$OMP PARALLEL 
943: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)   
944: !$OMP& REDUCTION(+:etemp,gradOMP) 
945:         etemp=0 
946:         TID=0 
947:         NTHREADS=1; 
948: !$      TID = OMP_GET_THREAD_NUM() 
949: !$      NTHREADS = OMP_GET_NUM_THREADS() 
950:         do i=1+TID,NATOMS,NTHREADS 
951: 926: 
952: !!!        do i=1,NATOMS927: !!        do i=1,NATOMS-1
953:            C1 = i928: !!                do j=i+1,NATOMS
954:                 Xgrid=int((X(i)-minX)/gridsize)+1 
955:                 Ygrid=int((Y(i)-minY)/gridsize)+1 
956:                 Zgrid=int((Z(i)-minZ)/gridsize)+1 
957:          
958:            do ii=XGRID-1,XGRID+1 
959:             do jj=YGRID-1,YGRID+1 
960:              do kk=ZGRID-1,ZGRID+1 
961:            if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and. 
962:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then 
963: 929: 
964:            do jjj=1,Ngrid(ii,jj,kk)930:            tempN=0
 931:         do i=1,maxgridX
 932:          do j=1,maxgridY
 933:           do k=1,maxgridZ
 934: 
 935:            do ii=i-1,i+1
 936:             do jj=j-1,j+1
 937:              do kk=k-1,k+1
 938:            if(i .ge. 1 .and. j .ge. 1 .and. k .ge. 1)then
 939:            if(i .eq. ii .and. j .eq. jj .and. k .eq. kk)then
 940: 
 941: 
 942:           do iii=1,Ngrid(i,j,k)-1
 943:            do jjj=iii+1,Ngrid(i,j,k)
 944: 
 945:            C1=grid(i,j,k,iii)
 946:            C2=grid(i,j,k,jjj)
 947:              if(NNCsigma(C1,C2) .gt. 0 .or. NNCsigma(C2,C1) .gt. 0)then
 948: ! only count each once
 949:                  tempN=tempN+1
 950:                  temparray(tempN,1)=C1
 951:                  temparray(tempN,2)=C2
 952:            endif
 953:            enddo
 954:           enddo
 955:  
 956:         else
965: 957: 
 958:           do iii=1,Ngrid(i,j,k)
 959:            do jjj=1,Ngrid(ii,jj,kk)
 960:            C1=grid(i,j,k,iii)
966:            C2=grid(ii,jj,kk,jjj)961:            C2=grid(ii,jj,kk,jjj)
967:              if(NNCinc(i,C2) .eq. 1)then 
968: !!! SSS 
969:         dx = X(C1) - X(C2) 
970:         dy = Y(C1) - Y(C2) 
971:         dz = Z(C1) - Z(C2) 
972:           r2 = dx**2 + dy**2 + dz**2 
973:           if(r2 .le. Rcut2)then 
974:              rm2 = 1/r2 
975:              rm14 = rm2**7 
976:                 etemp=etemp+STT*rm2**6+SC 
977:                 f_over_r=-STT*12.0*rm14 
978:                 RD1=sqrt(r2)-NCswitch 
979:                 if(r2 .gt. Rswitch2)then 
980:                         f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2) 
981:                         etemp=etemp+SA/3.0*RD1**3+SB/4.0*RD1**4 
982:                 elseif(r2 .lt. Rswitch2)then 
983:                 ! normal repulsive term 
984:                 else 
985:                 ! things should have fallen in one of the previous two... 
986:                 write(*,*) 'something went wrong with switching function' 
987:                 call abort 
988:                 endif 
989:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx 
990:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy 
991:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz 
992:  
993:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx 
994:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy 
995:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz 
996: 962: 
997:             endif963:              if(NNCsigma(C1,C2) .gt. 0 .or. NNCsigma(C2,C1) .gt. 0)then
 964:                  tempN=tempN+1
 965:                  temparray(tempN,1)=C1
 966:                  temparray(tempN,2)=C2
 967:            endif
998: 968: 
999: !!! END969:            enddo
 970:           enddo
1000: 971: 
1001:            endif 
1002:            enddo 
1003:   
1004:           endif972:           endif
1005: 973: 
 974:         endif
1006:             enddo975:             enddo
1007:            enddo976:            enddo
1008:           enddo977:           enddo
1009: 978: 
1010:           enddo979:           enddo
1011: !$OMP END PARALLEL 
1012:                 energy = energy  + etemp 
1013:         do C1=1,NATOMS 
1014:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2)  
1015:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1)  
1016:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)    
1017:         enddo 
1018:  
1019:       END 
1020:  
1021: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^ 
1022:  
1023: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
1024: ! FUNCTIONS for DH Calculations 
1025: !***************************************************** 
1026:  
1027:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r) 
1028:       USE KEY 
1029:       implicit NONE 
1030:       double precision kappa,r 
1031:       SBMDHV=exp(-kappa*r)/r 
1032:       END 
1033:  
1034:       DOUBLE PRECISION FUNCTION SBMDHVP(kappa,r) 
1035:       USE KEY 
1036:       implicit NONE 
1037:       double precision kappa,r 
1038:       SBMDHVP=-kappa*exp(-kappa*r)/r-exp(-kappa*r)/(r**2) 
1039:       END 
1040: !**************************************************** 
1041:  
1042: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
1043: !* SBMDHELEC computes the forces due to DH electrostatics        * 
1044: !********************************************************************** 
1045:  
1046:       subroutine SBMDHELEC(x,y,z,grad, energy, natoms,  
1047:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON) 
1048:  
1049:       USE KEY 
1050:       implicit NONE 
1051:       integer I, N, J, AN, NATOMS 
1052:  
1053:  
1054:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
1055:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC, 
1056:      Q SBMDHVP, SBMDHV 
1057:       integer C1, C2, ii,jj,kk,k,l,iii,jjj 
1058:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut  
1059:       DOUBLE PRECISION A,B,C, D, COEF2, COEF3 
1060:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS, 
1061:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS 
1062:  
1063:         INTEGER NNCinc(NATOMS,NATOMS) 
1064:         DOUBLE PRECISION dx,dy,dz 
1065:         integer tempN, alpha,SBMCHARGEON(NATOMS) 
1066:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS) 
1067:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut2 
1068:         ! Ngrid is the number of atoms in that grid point 
1069:         ! grid is the array of atoms in each grid 
1070:         integer Ngrid,grid,maxgrid,maxpergrid 
1071:         ! number of atoms per grid, max 
1072:         parameter (maxpergrid=50) 
1073:         ! dimensions of grid 
1074:         parameter (maxgrid=15) 
1075:         dimension Ngrid(maxgrid,maxgrid,maxgrid), 
1076:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid) 
1077:         integer MaxGridX,MaxGridY,MaxGridZ 
1078:         double precision gridsize,RD1 
1079:         double precision minX,minY,minZ,maxX,maxY,maxZ 
1080:         integer Xgrid,Ygrid,Zgrid,TID 
1081:         double precision Etemp,C2T 
1082:         diff=DHcut-DHswitch 
1083:         alpha=12 
1084:  
1085:  
1086:         GRIDSIZE=DHCUT*1.01 
1087:         DHcut2=DHcut**2 
1088:         DHswitch2=DHswitch**2 
1089:  
1090:         A=-PREFACTOR*SBMDHV(kappa,DHcut) 
1091:         B=-PREFACTOR*SBMDHVP(kappa,DHcut) 
1092:  
1093:         COEF3=(B-2*A/diff)/(diff**2)         
1094:         COEF2=(B-3*diff**2*COEF3)/(2*diff)         
1095: ! add DH interactions 
1096:  
1097: ! grid the system 
1098:         minX=10000000 
1099:         minY=10000000 
1100:         minZ=10000000 
1101:          
1102:         maxX=-10000000 
1103:         maxY=-10000000 
1104:         maxZ=-10000000 
1105:  
1106:         do ii=2,SBMCHARGEON(1) 
1107:            i=SBMCHARGEON(ii) 
1108:            gradOMP(i*3)=0 
1109:            gradOMP(i*3-1)=0 
1110:            gradOMP(i*3-2)=0 
1111:  
1112:            if(X(i) .gt. maxX)then 
1113:                 maxX=X(i) 
1114:            elseif(Y(i) .gt. maxY)then 
1115:                 maxY=Y(i) 
1116:            endif 
1117:  
1118:            if(Z(i) .gt. maxZ)then 
1119:                 maxZ=Z(i) 
1120:            elseif(X(i) .lt. minX)then 
1121:                 minX=X(i) 
1122:            endif 
1123:  
1124:            if(Y(i) .lt. minY)then 
1125:                 minY=Y(i) 
1126:            elseif(Z(i) .lt. minZ)then 
1127:                 minZ=Z(i) 
1128:            endif 
1129:         enddo 
1130:  
1131:         maxgridX=int((maxX-minX)/gridsize)+1 
1132:         maxgridY=int((maxY-minY)/gridsize)+1 
1133:         maxgridZ=int((maxZ-minZ)/gridsize)+1 
1134:  
1135:         if(maxgridX .ge. maxgrid .or.  
1136:      Q  maxgridY .ge. maxgrid .or. 
1137:      Q  maxgridZ .ge. maxgrid )then 
1138:         write(*,*) 'system got too big for grid searching...' 
1139:         call abort 
1140:         endif 
1141:  
1142:  
1143: !!  Add a second grid that only includes atoms that are charged 
1144:  
1145:         do i=1,maxgrid 
1146:          do j=1,maxgrid 
1147:           do k=1,maxgrid 
1148:                 Ngrid(i,j,k)=0 
1149:           enddo980:           enddo
1150:          enddo981:           enddo
1151:         enddo982: !        write(*,*) 'list prepared'
1152: 983:         !write(*,*) tempN, 'pairs'
1153:  
1154:         do ii=2,SBMCHARGEON(1) 
1155:                 i=SBMCHARGEON(ii) 
1156: 984: 
1157:                 Xgrid=int((X(i)-minX)/gridsize)+1985: !end of making list
1158:                 Ygrid=int((Y(i)-minY)/gridsize)+1986:            do i=1,tempN
1159:                 Zgrid=int((Z(i)-minZ)/gridsize)+1987:            
1160:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1988:            C1 = temparray(i,1)
1161:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then989:            C2 = temparray(i,2)
1162:                         write(*,*) 'too many atoms in a grid'990:            !ST=NNCsigma(C1,C2)
1163:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid991: 
1164:                         call abort992:         if(used(C1,C2) .neqv. .TRUE.)then
1165:                 endif993:         used(C1,C2) = .TRUE.
1166:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i994:         used(C2,C1) = .TRUE.
1167:         enddo 
1168: 995: 
 996:         dx = X(C1) - X(C2)
1169: 997: 
1170:            tempN=0998:          dy = Y(C1) - Y(C2)
1171: 999: 
 1000:         dz = Z(C1) - Z(C2)
1172: 1001: 
1173: !$OMP PARALLEL1002:           r2 = dx**2 + dy**2 + dz**2
1174: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)   
1175: !$OMP& REDUCTION(+:etemp,gradOMP) 
1176:         etemp=0 
1177:         etemp=01003:         etemp=0
1178:         TID=0 
1179:         NTHREADS=1; 
1180: !$      TID = OMP_GET_THREAD_NUM() 
1181: !$      NTHREADS = OMP_GET_NUM_THREADS() 
1182:         do iii=2+TID,SBMCHARGEON(1),NTHREADS 
1183: 1004: 
 1005:           if(r2 .le. Rcut2)then
 1006:              rm2 = 1/r2
 1007:              rm14 = rm2**7
1184: 1008: 
1185:            C1=SBMCHARGEON(iii)1009:                 etemp=STT*rm2**6+SC
1186:                 Xgrid=int((X(C1)-minX)/gridsize)+11010:                 f_over_r=-STT*12.0*rm14
1187:                 Ygrid=int((Y(C1)-minY)/gridsize)+1 
1188:                 Zgrid=int((Z(C1)-minZ)/gridsize)+1 
1189:          
1190:            do ii=XGRID-1,XGRID+1 
1191:             do jj=YGRID-1,YGRID+1 
1192:              do kk=ZGRID-1,ZGRID+1 
1193:            if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and. 
1194:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then 
1195: 1011: 
1196:            do jjj=1,Ngrid(ii,jj,kk)1012: !NCsigma(i), NNCeps(i)
 1013:                 RD1=sqrt(r2)-NCswitch
 1014:                 if(r2 .gt. Rswitch2)then
1197: 1015: 
1198:            C2=grid(ii,jj,kk,jjj)1016:                         f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
1199:            if(NNCinc(C1,C2) .eq. 1)then1017:                         etemp=etemp+SA/3.0*RD1**3+SB/4.0*RD1**4
1200: 1018: 
1201:         dx = X(C1) - X(C2)1019:                 elseif(r2 .lt. Rswitch2)then
1202:         dy = Y(C1) - Y(C2) 
1203:         dz = Z(C1) - Z(C2) 
1204:           r2 = dx**2 + dy**2 + dz**2 
1205:           if(r2 .le. DHcut2)then 
1206:            C2T=SBMCHARGE(C1)*SBMCHARGE(C2) 
1207:         ! add force evaluation now 
1208: 1020: 
1209:              rm2 = 1/r21021:                 ! normal repulsive term
1210:              r1=sqrt(r2) 
1211:                 etemp=etemp+PREFACTOR*C2T*SBMDHV(kappa,r1) 
1212: 1022: 
1213:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1) 
1214:                 if(r2 .gt. DHswitch2)then 
1215:                 RD1=r1-DHswitch 
1216:                         f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2)) 
1217:                         etemp=etemp+COEF2*RD1**2+COEF3*RD1**3 
1218:                 elseif(r2 .lt. DHswitch2)then 
1219:                 ! normal DH term 
1220:                 else1023:                 else
1221:                 ! things should have fallen in one of the previous two...1024:                 ! things should have fallen in one of the previous two...
1222:                 write(*,*) 'something went wrong with DH switching function'1025:                 write(*,*) 'something went wrong with switching function'
1223:                 call abort1026: !        call abort
1224:                 endif1027:                 endif
1225:                 f_over_r=f_over_r*sqrt(rm2) 
1226:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx 
1227:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy 
1228:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz 
1229:  
1230:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx 
1231:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy 
1232:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz 
1233: 1028: 
1234:             endif1029:                 energy = energy  + etemp
1235: 1030: 
1236:            endif 
1237:            enddo 
1238:   
1239:           endif 
1240: 1031: 
1241:             enddo1032: !                if(r2 .gt. Rswitch2*0.999 .and. r2 .lt. Rswitch2*1.001)then
1242:            enddo1033:                 !write(*,*) 'etemp',sqrt(r2),f_over_r,etemp
1243:           enddo1034: !                endif
1244: 1035: 
1245:           enddo 
1246: !$OMP END PARALLEL 
1247:                 energy = energy  + etemp 
1248:         do i=2,SBMCHARGEON(1) 
1249:            C1=SBMCHARGEON(i) 
1250:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2)  
1251:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1)  
1252:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)    
1253:         enddo 
1254: 1036: 
1255:       END1037: ! f_over_r is the force over the magnitude of r so there is no need to resolve
 1038: ! the dx, dy and dz into unit vectors
1256: 1039: 
1257: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^1040: ! now add the acceleration 
 1041:               grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx
 1042:               grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy
 1043:               grad(C1*3)   = grad(C1*3)   + f_over_r * dz
 1044: 
 1045:                grad(C2*3-2) =  grad(C2*3-2) - f_over_r * dx
 1046:               grad(C2*3-1) =  grad(C2*3-1) - f_over_r * dy
 1047:               grad(C2*3)   =  grad(C2*3)  -  f_over_r * dz
 1048:             endif
1258: 1049: 
1259: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1050:             endif
1260: !* SBMPR computes the forces due to position restraints               * 
1261: !********************************************************************** 
1262: 1051: 
1263:       subroutine SBMPR(x,y,z,grad, energy, natoms, 1052:             enddo
1264:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 
1265: 1053: 
1266:       USE KEY1054:            do i=1,tempN
1267:       implicit NONE1055:            
1268:       integer I, J, NATOMS1056:            C1 = temparray(i,1)
1269:       INTEGER SBMPRN, SBMPRI(NATOMS)1057:            C2 = temparray(i,2)
1270:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)1058:            !ST=NNCsigma(C1,C2)
1271:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),1059: 
1272:      Q grad(3*NATOMS), energy1060:                 used(C1,C2) = .FALSE.
1273:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K61061:                 used(C2,C1) = .FALSE.
1274:               do I=1,SBMPRN1062:            enddo
1275:            J=SBMPRI(I)1063: 
1276:            DX=X(J)-SBMPRX(I,1)1064: 
1277:            DY=Y(J)-SBMPRX(I,2)1065: !!           end do
1278:            DZ=Z(J)-SBMPRX(I,3)1066: !!        enddo
1279:            K1=SBMPRK(I,1) 
1280:            K2=SBMPRK(I,2) 
1281:            K3=SBMPRK(I,3) 
1282:            K4=SBMPRK(I,4) 
1283:            K5=SBMPRK(I,5) 
1284:            K6=SBMPRK(I,6) 
1285:  
1286:  
1287:            ENERGY = ENERGY +  
1288:      Q     0.5*K1*DX**2 +  
1289:      Q     0.5*K2*DY**2 + 
1290:      Q     0.5*K3*DZ**2 + 
1291:      Q     K4*DX*DY + 
1292:      Q     K5*DX*DZ + 
1293:      Q     K6*DY*DZ 
1294:  
1295:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ  
1296:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ 
1297:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY 
1298:         enddo 
1299: 1067: 
1300:       END1068:       END
1301: 1069: 
1302: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^1070: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
 1071: 
1303: 1072: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0