hdiff output

r30668/SBM.f 2016-07-06 15:38:16.405399242 +0100 r30667/SBM.f 2016-07-06 15:38:16.741403790 +0100
  1:       SUBROUTINE SBM(qo,NATOMS,grad,energy,GTEST,STEST)  1:       SUBROUTINE SBM(qo,NATOMS,grad,energy,GTEST,STEST)
  2:       USE KEY  2:       USE KEY
  3:       implicit NONE  3:       implicit NONE
  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=25000) 11:         parameter(NSBMMAX=5000)
 12:       INTEGER MAXSEP,MAXSEP1 12:       DOUBLE PRECISION  Rb(NSBMMAX), bK(NSBMMAX), ANTC(NSBMMAX),
 13:       PARAMETER (MAXSEP1=30) 13:      Q Tk(NSBMMAX), PK(NSBMMAX), PHISBM(NSBMMAX), Sigma(NSBMMAX*5),
 14:       LOGICAL NNCinc(NSBMMAX,MAXSEP1) 14:      Q EpsC(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
 15:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 15:       INTEGER  Ib1(NSBMMAX), Ib2(NSBMMAX), IT(NSBMMAX), JT(NSBMMAX),
 16:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF1(NSBMMAX*5), 16:      Q KT(NSBMMAX),IP(NSBMMAX), JP(NSBMMAX), KP(NSBMMAX),
 17:      Q CONCOEF2(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 17:      Q LP(NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5),
 18:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2), 18:      Q PHITYPE(NSBMMAX),NBA,NTA,NPA,NC,NNCinc(NSBMMAX,NSBMMAX),SBMPRN, SBMPRI(NSBMMAX)
 19:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3), 
 20:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5), 
 21:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX) 
 22:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3) 19:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3)
 23:         integer CONTACTTYPE,SBMCHARGEON(NSBMMAX) 20:         integer CONTACTTYPE,SBMCHARGEON(NSBMMAX)
 24:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 21:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 25:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 22:      Q sigma, epsC,NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
 26:      Q SBMPRK,SBMPRX 23:      Q SBMPRK,SBMPRX
 27:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP, 24:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP,
 28:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC, 25:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC,NNCinc,
 29:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,MAXSEP 26:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI
 30:         COMMON /logical/ NNCinc 
 31:         MAXSEP=MAXSEP1 
 32:         if(NATOMS.gt. NSBMMAX)then 27:         if(NATOMS.gt. NSBMMAX)then
 33:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 28:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 34:         STOP 29:         STOP
 35:         endif 30:         endif
 36:  31: 
 37:        if(.NOT.CALLED)then 32:        if(.NOT.CALLED)then
 38:         write(*,*) 33:         write(*,*)
 39:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 34:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 40:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 35:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 41:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 36:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 42:         write(*,*) 37:         write(*,*)
 43:          38:         
 44:  39: 
 45:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP,  40:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP, 
 46:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2,  41:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, Sigma, 
 47:      Q NNCinc,MAXSEP,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT, 42:      Q EpsC, NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,
 48:      Q NSBMMAX,CONTACTTYPE, 43:      Q NSBMMAX,CONTACTTYPE,
 49:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 44:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 50:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or. 45:         if(NBA .gt. NSBMMAX .or. NTA .gt. NSBMMAX .or.
 51:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then 46:      Q  NPA .gt. NSBMMAX .or. NC .gt. 5*NSBMMAX)then
 52:         write(*,*) 'increase array size' 47:         write(*,*) 'increase array size'
 53:         write(*,*) NBA, NTA,NPA,NC  
 54:         stop 48:         stop
 55:         endif  49:         endif 
  50: 
 56:         CALLED=.TRUE. 51:         CALLED=.TRUE.
 57:         endIF 52:         endIF
 58: ! call the energy routine 53: ! call the energy routine
 59:  54: 
 60:       call calc_energy_SBM(qo,natoms,GRAD,energy,Ib1,Ib2, 55:       call calc_energy_SBM(qo,natoms, GRAD, energy, Ib1, Ib2,
 61:      Q Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK, 56:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,
 62:      Q PHITYPE,PHISBM,IC,JC,CONCOEF1,CONCOEF2,  57:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC, 
 63:      Q  NNCinc,MAXSEP,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON, 58:      Q  NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,
 64:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 59:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
 65:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 60:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 66:       IF (STEST) THEN 61:       IF (STEST) THEN
 67:          PRINT '(A)','ERROR - second derivatives not available' 62:          PRINT '(A)','ERROR - second derivatives not available'
 68:          STOP 63:          STOP
 69:       ENDIF 64:       ENDIF
 70:       return 65:       return
 71:       end 66:       end
 72:  67: 
 73:  68: 
 74: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 69: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 75: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 70: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 76: !* startt then the velocities are assigned, otherwise, they are read   * 71: !* startt then the velocities are assigned, otherwise, they are read   *
 77: !* by selecting 2, or generated by selecting 3                         * 72: !* by selecting 2, or generated by selecting 3                         *
 78: !*********************************************************************** 73: !***********************************************************************
 79:  74: 
 80:       subroutine SBMinit(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk, 75:       subroutine SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk,
 81:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2,  76:      Q  IP, JP, KP, LP, PK,PHITYPE,PHISBM,IC, JC, Sigma, 
 82:      Q NNCinc,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch, 77:      Q EpsC, NNCinc, NBA, 
 83:      Q NCcut,STT,NSBMMAX,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 78:      Q NTA, NPA, NC, SBMCHARGE, SBMCHARGEON,NCswitch,NCcut,STT,NSBMMAX,
  79:      Q CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
 84:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 80:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 85:       USE KEY 81:       USE KEY
 86:       USE COMMONS, only: ATMASS 82:       USE COMMONS, only: ATMASS
 87:       implicit NONE 83:       implicit NONE
 88:  84: 
 89:         integer i,j,MaxCon,NNCmax,NATOMS,storage, dummy,  ANr,  85:         integer i,j,MaxCon,NNCmax,NATOMS,storage, dummy,  ANr, IB11,
 90:      Q  Ib22, Ib21,IT1, JT1, KT1, IT2, JT2, KT2, IP1, JP1, 86:      Q IB12, Ib22, Ib21,IT1, JT1, KT1, IT2, JT2, KT2, IP1, JP1,
 91:      Q KP1, LP1, IP2, JP2, KP2, 87:      Q KP1, LP1, IP2, JP2, KP2,
 92:      Q LP2, nBA1, nTA1, nPA1, nBA2, nTA2, nPA2,  ind1, ind2, ANt, 88:      Q LP2, nBA1, nTA1, nPA1, nBA2, nTA2, nPA2,  ind1, ind2, ANt,
 93:      Q  MDT1, MDT2, cl1, cl2,tempi,NSBMMAX,MAXSEP,MAXSEPTEMP 89:      Q  MDT1, MDT2, cl1, cl2,tempi,NSBMMAX
 94:  90: 
 95:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX),  91:       DOUBLE PRECISION  Rb(NSBMMAX), bK(NSBMMAX), ANTC(NSBMMAX), 
 96:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), 92:      Q Tk(NSBMMAX), PK(NSBMMAX), PHISBM(NSBMMAX),
 97:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS), 93:      Q Sigma(NSBMMAX*5), EpsC(NSBMMAX*5),  
 98:      Q CONCOEF1(5*NSBMMAX),CONCOEF2(5*NSBMMAX)  94:      Q NCswitch,NCcut,STT,SBMCHARGE(NATOMS)
 99:       INTEGER  Ib1(2*NSBMMAX), Ib2(2*NSBMMAX), IT(2*NSBMMAX), JT(2*NSBMMAX),  95:       INTEGER  Ib1(NSBMMAX), Ib2(NSBMMAX), IT(NSBMMAX), JT(NSBMMAX), 
100:      Q KT(2*NSBMMAX),IP(3*NSBMMAX), JP(3*NSBMMAX), KP(3*NSBMMAX), 96:      Q KT(NSBMMAX),IP(NSBMMAX), JP(NSBMMAX), KP(NSBMMAX),
101:      Q LP(3*NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5),  97:      Q LP(NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 
102:      Q  NBA, NTA, NPA, NC,SBMCHARGEON(NATOMS) 98:      Q  NBA, NTA, NPA, NC,SBMCHARGEON(NATOMS)
103:       INTEGER PHITYPE(3*NSBMMAX) 99:       INTEGER PHITYPE(NSBMMAX),NNCinc(NATOMS,NATOMS)
 100:        DOUBLE PRECISION  pinitmax, TK1, TK2,  APTtemp, msT,
 101:      Q SigmaT1, SigmaT2, epstemp
104:       integer AA,BB,ANTEMP102:       integer AA,BB,ANTEMP
105:       LOGICAL NNCINC(NATOMS,MAXSEP) 
106:       DOUBLE PRECISION dx,dy,dz103:       DOUBLE PRECISION dx,dy,dz
107:        double precision PI104:        double precision PI
108:         DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut105:         DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
 106:         logical TEMPARRAY
 107:         dimension TEMPARRAY(NATOMS,NATOMS)
109:         integer CONTACTTYPE108:         integer CONTACTTYPE
110:         character TMPCHAR109:         character TMPCHAR
111:         integer TMPINT,NUMCHARGES,nexc,I1,I2110:         integer TMPINT,NUMCHARGES,nexc,I1,I2
112:         double precision TMPREAL,concentration111:         double precision TMPREAL,concentration
113: 112: 
114:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)113:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)
115:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)114:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
116: 115: 
117: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 116: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
118: !$OMP PARALLEL117: !$OMP PARALLEL
119: !$    NTHREADS = OMP_GET_NUM_THREADS()118: !$    NTHREADS = OMP_GET_NUM_THREADS()
120: !$      TID = OMP_GET_THREAD_NUM()119: !$      TID = OMP_GET_THREAD_NUM()
121: !$    if(TID .eq. 0)then 120: !$    if(TID .eq. 0)then 
122: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS121: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS
123: !$    endif122: !$    endif
124: !$OMP END PARALLEL123: !$OMP END PARALLEL
125:       pi = 3.14159265358979323846264338327950288419716939937510124:       pi = 3.14159265358979323846264338327950288419716939937510
126:       MAXSEPTEMP=0125: 
127:       NNCmax = NATOMS*NATOMS126:       NNCmax = NATOMS*NATOMS
128:       MaxCon=NATOMS*5127:       MaxCon=NATOMS*5
129: 128: 
130:       do i=1,NATOMS129:       do i=1,NATOMS
131:         do j=1,30130:         do j=1,NATOMS
132:           NNCINC(i,j)=.TRUE.131:           TEMPARRAY(i,j)=.TRUE.
133:         enddo132:         enddo
134:       enddo133:       enddo
135: 134: 
136: 135: 
137: 136: 
138: ! These lines read in the parameters.137: ! These lines read in the parameters.
139:         open(30, file='SBM.INP', status='old', access='sequential')138:         open(30, file='SBM.INP', status='old', access='sequential')
140:         write(*,*) 139:         write(*,*) 
141:         write(*,*) 'Reading the forcefield from SBM.INP'140:         write(*,*) 'Reading the forcefield from SBM.INP'
142:         write(*,*) 141:         write(*,*) 
168:         read(30,*) 167:         read(30,*) 
169:         read(30,*) NC168:         read(30,*) NC
170:         write(*,*) NC, 'contacts'        169:         write(*,*) NC, 'contacts'        
171:         read(30,*) CONTACTTYPE170:         read(30,*) CONTACTTYPE
172:           if(NC .gt. MaxCon)then171:           if(NC .gt. MaxCon)then
173:              write(*,*) 'too many contacts'172:              write(*,*) 'too many contacts'
174:              STOP173:              STOP
175:           endif174:           endif
176: 175: 
177:         do i=1, NC176:         do i=1, NC
178:           read(30, *) IC(i), JC(i), Sigma, EpsC177:           read(30, *) IC(i), JC(i), Sigma(i), EpsC(i)
179: 178:           TEMPARRAY(IC(i), JC(i))=.FALSE.
180:           if(CONTACTTYPE .eq. 1)then179:           TEMPARRAY(JC(i), IC(i))=.FALSE.
181:             ConCoef1(I)=EPSC*SIGMA**12-REPS*RSIG**12180:           Sigma(i)=Sigma(i)**2
182:             ConCoef2(I)=2*EPSC*SIGMA**6 
183:           elseif(CONTACTTYPE .eq. 2)then 
184:             ConCoef1(I)=5*EPSC*SIGMA**12-REPS*RSIG**12 
185:             ConCoef2(I)=6*EPSC*SIGMA**10 
186:           else 
187:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.' 
188:           endif 
189: !          AA=min(IC(I),JC(I)) 
190: !          BB=max(IC(I),JC(I)) 
191: !          BB=BB-AA 
192: !          if(BB .gt. MAXSEP)then 
193: !                write(*,*) 'ERROR: exclusions between atoms too far apart. Sep=',BB 
194: !          endif 
195: !          NNCINC(AA, BB)=.FALSE. 
196: ! NEED TO SUBTRACT OUT THE EXCLUDED VOLUME EFFECT,CHANGE TO 6-12, or 
197: ! 10-12 coefficients 
198:         end do181:         end do
199: 182: 
200:           read(30,*)183:           read(30,*)
201:           read(30,*) nBA184:           read(30,*) nBA
202:           write(*,*) nBA, 'bonds'185:           write(*,*) nBA, 'bonds'
203:         do i=1, nBA186:         do i=1, nBA
204:           read(30,*) Ib1(i), Ib2(i),Rb(i), bK(i)187:           read(30,*) Ib1(i), Ib2(i),Rb(i), bK(i)
205:           AA=min(IB1(I),IB2(I))188:           TEMPARRAY(Ib1(i), Ib2(i))=.FALSE.
206:           BB=max(IB1(I),IB2(I))189:           TEMPARRAY(Ib2(i), Ib1(i))=.FALSE.
207:           BB=BB-AA 
208:           if(BB .gt. MAXSEPTEMP)then 
209:             MAXSEPTEMP=BB 
210:           endif 
211:           if(BB .gt. MAXSEP)then 
212:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
213:           endif 
214:           NNCINC(AA, BB)=.FALSE. 
215:  
216:         end do190:         end do
217: 191: 
218:           read(30,*)192:           read(30,*)
219:           read(30,*) nTA193:           read(30,*) nTA
220:           write(*,*) nTA, 'angles'194:           write(*,*) nTA, 'angles'
221:         do i=1, nTA195:         do i=1, nTA
222:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)196:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)
223: 197:           TEMPARRAY(IT(i), JT(i))=.FALSE.
224:           AA=min(IT(I),JT(I))198:           TEMPARRAY(JT(i), IT(i))=.FALSE.
225:           BB=max(IT(I),JT(I))199:           TEMPARRAY(IT(i), KT(i))=.FALSE.
226:           BB=BB-AA200:           TEMPARRAY(KT(i), IT(i))=.FALSE.
227:           if(BB .gt. MAXSEP)then201:           TEMPARRAY(JT(i), KT(i))=.FALSE.
228:                 write(*,*) 'ERROR: exclusions between atoms too far apart',AA,BB202:           TEMPARRAY(KT(i), JT(i))=.FALSE.
229:           endif 
230:           NNCINC(AA, BB)=.FALSE. 
231:           AA=min(IT(I),KT(I)) 
232:           BB=max(IT(I),KT(I)) 
233:           BB=BB-AA 
234:           if(BB .gt. MAXSEPTEMP)then 
235:            MAXSEPTEMP=BB 
236:           endif 
237:           if(BB .gt. MAXSEP)then 
238:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
239:           endif 
240:           NNCINC(AA, BB)=.FALSE. 
241:           AA=min(KT(I),JT(I)) 
242:           BB=max(KT(I),JT(I)) 
243:           BB=BB-AA 
244:           if(BB .gt. MAXSEPTEMP)then 
245:            MAXSEPTEMP=BB 
246:           endif 
247:  
248:           if(BB .gt. MAXSEP)then 
249:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
250:           endif 
251:           NNCINC(AA, BB)=.FALSE. 
252:         enddo203:         enddo
253: 204: 
254:           read(30,*) 205:           read(30,*) 
255:           read(30,*) nPA206:           read(30,*) nPA
256:           write(*,*) nPA, 'dihedrals'207:           write(*,*) nPA, 'dihedrals'
257: 208: 
258: ! this reads in the dihedral angles and calculates the cosines and sines209: ! this reads in the dihedral angles and calculates the cosines and sines
259: ! in order to make the force and energy calculations easier, later.210: ! in order to make the force and energy calculations easier, later.
260:         do i=1, npA211:         do i=1, npA
261:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)212:            read(30,*) IP(i),JP(i),KP(i),LP(i),
262:           AA=min(IP(I),JP(I))213:      Q  PHITYPE(i),PHISBM(i),PK(i)
263:           BB=max(IP(I),JP(I))214:           TEMPARRAY(IP(i), JP(i))=.FALSE.
264:           BB=BB-AA215:           TEMPARRAY(JP(i), IP(i))=.FALSE.
265:           if(BB .gt. MAXSEPTEMP)then216:           TEMPARRAY(IP(i), KP(i))=.FALSE.
266:            MAXSEPTEMP=BB217:           TEMPARRAY(KP(i), IP(i))=.FALSE.
267:           endif218:           TEMPARRAY(IP(i), LP(i))=.FALSE.
268: 219:           TEMPARRAY(LP(i), IP(i))=.FALSE.
269:           if(BB .gt. MAXSEP)then220:           TEMPARRAY(JP(i), KP(i))=.FALSE.
270:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB221:           TEMPARRAY(KP(i), JP(i))=.FALSE.
271:           endif222:           TEMPARRAY(JP(i), LP(i))=.FALSE.
272:           NNCINC(AA, BB)=.FALSE.223:           TEMPARRAY(LP(i), JP(i))=.FALSE.
273:           AA=min(IP(I),KP(I))224:           TEMPARRAY(KP(i), LP(i))=.FALSE.
274:           BB=max(IP(I),KP(I))225:           TEMPARRAY(LP(i), KP(i))=.FALSE.
275:           BB=BB-AA226:         
276:           if(BB .gt. MAXSEPTEMP)then 
277:            MAXSEPTEMP=BB 
278:           endif 
279:  
280:           if(BB .gt. MAXSEP)then 
281:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
282:           endif 
283:           NNCINC(AA, BB)=.FALSE. 
284:           AA=min(IP(I),LP(I)) 
285:           BB=max(IP(I),LP(I)) 
286:           BB=BB-AA 
287:           if(BB .gt. MAXSEPTEMP)then 
288:            MAXSEPTEMP=BB 
289:           endif 
290:  
291:           if(BB .gt. MAXSEP)then 
292:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
293:           endif 
294:           NNCINC(AA, BB)=.FALSE. 
295:           AA=min(JP(I),KP(I)) 
296:           BB=max(JP(I),KP(I)) 
297:           BB=BB-AA 
298:           if(BB .gt. MAXSEPTEMP)then 
299:            MAXSEPTEMP=BB 
300:           endif 
301:  
302:           if(BB .gt. MAXSEP)then 
303:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
304:           endif 
305:           NNCINC(AA, BB)=.FALSE. 
306:           AA=min(JP(I),LP(I)) 
307:           BB=max(JP(I),LP(I)) 
308:           BB=BB-AA 
309:           if(BB .gt. MAXSEPTEMP)then 
310:            MAXSEPTEMP=BB 
311:           endif 
312:  
313:           if(BB .gt. MAXSEP)then 
314:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
315:           endif 
316:           NNCINC(AA, BB)=.FALSE. 
317:           AA=min(KP(I),LP(I)) 
318:           BB=max(KP(I),LP(I)) 
319:           BB=BB-AA 
320:           if(BB .gt. MAXSEPTEMP)then 
321:            MAXSEPTEMP=BB 
322:           endif 
323:  
324:           if(BB .gt. MAXSEP)then 
325:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
326:           endif 
327:           NNCINC(AA, BB)=.FALSE. 
328:         END DO227:         END DO
329: 228: 
330:           read(30,*)229:           read(30,*)
331:           read(30,*) nexc230:           read(30,*) nexc
332:           write(*,*) nexc, 'exculusions'231:           write(*,*) nexc, 'exculusions'
333: 232: 
334:         do i=1, nexc233:         do i=1, nexc
335:           read(30,*) I1, I2234:           read(30,*) I1, I2
336:           AA=min(I1,I2)235:           TEMPARRAY(I1, I2)=.FALSE.
337:           BB=max(I1,I2)236:           TEMPARRAY(I2, I1)=.FALSE.
338:           BB=BB-AA 
339:           if(BB .gt. MAXSEPTEMP)then 
340:            MAXSEPTEMP=BB 
341:           endif 
342:  
343:           if(BB .gt. MAXSEP)then 
344:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB 
345:           endif 
346:           NNCINC(AA, BB)=.FALSE. 
347:         end do237:         end do
348: 238: 
349: 239: 
350:         STT=reps*rsig**12240:         do i=1,NATOMS-1
 241:                 do j=i+1,NATOMS
 242:                 if(TEMPARRAY(i,j) .eqv. .TRUE.)then
 243: ! this simplifies calculations later
 244:          NNCinc(i,j) = 1
 245:          NNCinc(j,i) = 0
 246:                 else
 247:          NNCinc(i,j) = 0
 248:          NNCinc(j,i) = 0
 249: 
 250:                 endif
 251:                 enddo
 252:         end do
 253: 
 254:                 STT=reps*rsig**12
351: 255: 
352:         do I=1,NATOMS256:         do I=1,NATOMS
353:           TARR(I)=0257:           TARR(I)=0
354:         enddo258:         enddo
355: 259: 
356: ! read in position restraints260: ! read in position restraints
357:        read(30,*) 261:        read(30,*) 
358:        read(30,*) SBMPRN262:        read(30,*) SBMPRN
359:        write(*,*) SBMPRN, 'position restraints'263:        write(*,*) SBMPRN, 'position restraints'
360:        do I=1,SBMPRN264:        do I=1,SBMPRN
363:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)267:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)
364:             if(TARR(SBMPRI(I)) .ne. 0)then268:             if(TARR(SBMPRI(I)) .ne. 0)then
365:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)269:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
366: !                call abort270: !                call abort
367:             endif271:             endif
368:             TARR(SBMPRI(I))=1272:             TARR(SBMPRI(I))=1
369:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  273:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  
370:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'274:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'
371:         endif275:         endif
372:        enddo 276:        enddo 
373:         MAXSEP=MAXSEPTEMP277: 
374:        close(30)278:        close(30)
375:        end279:        end
376: 280: 
377: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^281: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
378: 282: 
379: 283: 
380: C284: C
381: C Calculate the Forces and energies285: C Calculate the Forces and energies
382: C286: C
383:       subroutine calc_energy_SBM(qo,natoms,GRAD, energy,Ib1, Ib2,287:       subroutine calc_energy_SBM(qo,natoms,GRAD, energy,Ib1, Ib2,
384:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,288:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,
385:      Q PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2, 289:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC, 
386:      Q NNCinc,MAXSEP,NBA, NTA, NPA, NC,SBMCHARGE, SBMCHARGEON,NCswitch,290:      Q NNCinc,NBA, NTA, NPA, NC,SBMCHARGE, SBMCHARGEON,NCswitch,
387:      Q STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,291:      Q STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
388:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)292:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
389: 293: 
390:       INTEGER I, J, NATOMS,NBA, NTA, NPA, NC,MAXSEP294:       INTEGER I, J, NATOMS,NBA, NTA, NPA, NC
391: 295: 
392:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY296:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
393:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)297:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
394: 298: 
395:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 299:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
396:      Q PHISBM(NPA),CONCOEF1(NC),CONCOEF2(NC),NCswitch,NCcut,STT, 300:      Q PHISBM(NPA), Sigma(NC),NCswitch,NCcut,STT, 
397:      Q SBMCHARGE(NATOMS)301:      Q EpsC(NC),SBMCHARGE(NATOMS)
398:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 302:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 
399:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),SBMCHARGEON(NATOMS),303:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),SBMCHARGEON(NATOMS),
400:      Q PHITYPE(NPA),CONTACTTYPE304:      Q PHITYPE(NPA),CONTACTTYPE,NNCinc(NATOMS,NATOMS)
401:         LOGICAL NNCinc(NATOMS,MAXSEP) 
402:       DOUBLE PRECISION dx,dy,dz305:       DOUBLE PRECISION dx,dy,dz
403:       INTEGER SBMPRN, SBMPRI(NATOMS)306:       INTEGER SBMPRN, SBMPRI(NATOMS)
404:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)307:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
405:       do i = 1, natoms308:       do i = 1, natoms
406:         j = (i-1)*3309:         j = (i-1)*3
407:         x(i) = qo(j+1)310:         x(i) = qo(j+1)
408:         y(i) = qo(j+2)311:         y(i) = qo(j+2)
409:         z(i) = qo(j+3)312:         z(i) = qo(j+3)
410:         grad(j+1) = 0.0313:         grad(j+1) = 0.0
411:         grad(j+2) = 0.0314:         grad(j+2) = 0.0
412:         grad(j+3) = 0.0315:         grad(j+3) = 0.0
413:       enddo316:       enddo
414: 317: 
415:       energy = 0.0318:       energy = 0.0
416:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)319:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
417:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)320:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
418:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,321:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
419:      Q PHITYPE,PHISBM,NPA)322:      Q PHITYPE,PHISBM,NPA)
420:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC, 323:       call SBMContacts(x,y,z,grad, energy, natoms, IC, JC, 
421:      Q CONCOEF1,CONCOEF2,NC,CONTACTTYPE)324:      Q Sigma, EpsC, NC,CONTACTTYPE)
422:       call SBMNonContacts(x,y,z,grad, energy, natoms, 325:       call SBMNonContacts(x,y,z,grad, energy, natoms, 
423:      Q NNCinc,MAXSEP,NCswitch,NCcut,STT)326:      Q NNCinc,NCswitch,NCcut,STT)
424:       call SBMDHELEC(x,y,z,grad, energy, natoms, 327:       call SBMDHELEC(x,y,z,grad, energy, natoms, 
425:      Q NNCinc,MAXSEP,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)328:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
426:       call SBMPR(x,y,z,grad, energy, natoms, 329:       call SBMPR(x,y,z,grad, energy, natoms, 
427:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)330:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
428: 331: 
429:       end332:       end
430: 333: 
431: 334: 
432: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<335: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
433: !* SBMBonds  computes the hookean force and energy between chosen atoms *336: !* SBMBonds  computes the hookean force and energy between chosen atoms *
434: !***********************************************************************337: !***********************************************************************
435: 338: 
719: 622: 
720: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^623: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^
721: 624: 
722: 625: 
723: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<626: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
724: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *627: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *
725: !* 10-12 or 6-12 potential                                              *628: !* 10-12 or 6-12 potential                                              *
726: !***********************************************************************629: !***********************************************************************
727: 630: 
728:       subroutine SBMcontacts(x,y,z,grad,energy,631:       subroutine SBMcontacts(x,y,z,grad,energy,
729:      Q NATOMS,IC,JC,CONCOEF1,CONCOEF2,NC,CONTACTTYPE)632:      Q NATOMS,IC,JC,Sigma,EpsC,NC,CONTACTTYPE)
730:       USE KEY633:       USE KEY
731:       implicit NONE634:       implicit NONE
732:       integer I, NATOMS,NC635:       integer I, NATOMS,NC
733: 636: 
734:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 637:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
735:      Q grad(3*NATOMS), energy638:      Q grad(3*NATOMS), energy
736:       DOUBLE PRECISION dx,dy,dz639:       DOUBLE PRECISION dx,dy,dz
737: 640: 
738:       integer C1, C2 641:       integer C1, C2 
739:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 642:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 
740: 643: 
741:         DOUBLE PRECISION CONCOEF1(NC), CONCOEF2(NC)644:         DOUBLE PRECISION Sigma(NC), EpsC(NC)
742:         INTEGER IC(NC), JC(NC),CONTACTTYPE645:         INTEGER IC(NC), JC(NC),CONTACTTYPE
743: 646: 
744: ! type 1 is 6-12647: ! type 1 is 6-12
745: ! type 2 is 12-12648: ! type 2 is 12-12
746: ! implement type 3 as gaussian649: ! implement type 3 as gaussian
747: ! implement type 4 as dual-basin gaussian650: ! implement type 4 as dual-basin gaussian
748: ! question: Should type be part of the interaction, or should they be organized651: ! question: Should type be part of the interaction, or should they be organized
749: ! into separate lists?652: ! into separate lists?
750: 653: 
751: ! type 1 is 6-12 interaction654: ! type 1 is 6-12 interaction
753: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R)REDUCTION(+:ENERGY,grad)656: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R)REDUCTION(+:ENERGY,grad)
754: !$OMP DO657: !$OMP DO
755:           do i=1, NC658:           do i=1, NC
756:             C1 = IC(i)659:             C1 = IC(i)
757:             C2 = JC(i)660:             C2 = JC(i)
758:             dx = X(C1) - X(C2)661:             dx = X(C1) - X(C2)
759:             dy = Y(C1) - Y(C2)662:             dy = Y(C1) - Y(C2)
760:             dz = Z(C1) - Z(C2)663:             dz = Z(C1) - Z(C2)
761:             r2 = dx**2 + dy**2 + dz**2664:             r2 = dx**2 + dy**2 + dz**2
762:             rm2 = 1.0/r2665:             rm2 = 1.0/r2
763:           !  rm2 = rm2*sigma(i)666:             rm2 = rm2*sigma(i)
764:             rm6 = rm2**3667:             rm6 = rm2**3
765:             !ENERGY = ENERGY + epsC(i)*rm6*(rm6-2.0)668:             ENERGY = ENERGY + epsC(i)*rm6*(rm6-2.0)
766:             ENERGY = ENERGY + rm6*(CONCOEF1(I)*rm6-CONCOEF2(I))669:             f_over_r = -epsC(i)*12.0*rm6*(rm6-1.0)/r2
767:             f_over_r = -12.0*rm6*(CONCOEF1(I)*rm6-CONCOEF2(I)/2.0)*rm2 
768:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx670:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
769:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy671:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
770:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz672:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz
771:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx673:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
772:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy674:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
773:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz675:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz
774:           enddo676:           enddo
775: !$OMP END DO677: !$OMP END DO
776: !$OMP END PARALLEL678: !$OMP END PARALLEL
777: ! type 2 is 10-12 interaction679: ! type 2 is 10-12 interaction
779: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R)REDUCTION(+:ENERGY,grad)681: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R)REDUCTION(+:ENERGY,grad)
780: !$OMP DO682: !$OMP DO
781:           do i=1, NC683:           do i=1, NC
782:              C1 = IC(i)684:              C1 = IC(i)
783:              C2 = JC(i)685:              C2 = JC(i)
784:              dx = X(C1) - X(C2)686:              dx = X(C1) - X(C2)
785:              dy = Y(C1) - Y(C2)687:              dy = Y(C1) - Y(C2)
786:              dz = Z(C1) - Z(C2)688:              dz = Z(C1) - Z(C2)
787:              r2 = dx**2 + dy**2 + dz**2689:              r2 = dx**2 + dy**2 + dz**2
788:              rm2 = 1.0/r2690:              rm2 = 1.0/r2
789: !             rm2 = rm2*sigma(i)691:              rm2 = rm2*sigma(i)
790:              rm10 = rm2**5692:              rm10 = rm2**5
791:              !ENERGY = ENERGY + epsC(i)*rm10*(5*rm2-6.0)693:              ENERGY = ENERGY + epsC(i)*rm10*(5*rm2-6.0)
792:              ENERGY = ENERGY + rm10*(CONCOEF1(I)*rm2-CONCOEF2(I))694:              f_over_r = -epsc(i)*60.0*rm10*(rm2-1.0)/r2
793:              f_over_r = -60.0*rm10*(CONCOEF1(I)*rm2/5.0-CONCOEF2(I)/6.0)*rm2 
794:              grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx695:              grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
795:              grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy696:              grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
796:              grad(3*C1)   = grad(3*C1)   + f_over_r * dz697:              grad(3*C1)   = grad(3*C1)   + f_over_r * dz
797:              grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx698:              grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
798:              grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy699:              grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
799:              grad(3*C2)   = grad(3*C2)   - f_over_r * dz700:              grad(3*C2)   = grad(3*C2)   - f_over_r * dz
800:            enddo701:            enddo
801: !$OMP END DO702: !$OMP END DO
802: !$OMP END PARALLEL703: !$OMP END PARALLEL
803: 704: 
810:       end711:       end
811: 712: 
812: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^713: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
813: 714: 
814: 715: 
815: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<716: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
816: !* SBMNonContacts computes the forces due to non native contacts       *717: !* SBMNonContacts computes the forces due to non native contacts       *
817: !**********************************************************************718: !**********************************************************************
818: 719: 
819:       subroutine SBMnoncontacts(x,y,z,grad, energy, 720:       subroutine SBMnoncontacts(x,y,z,grad, energy, 
820:      Q NATOMS,NNCinc,MAXSEP,NCswitch,NCcut,STT)721:      Q NATOMS,NNCinc,NCswitch,NCcut,STT)
821:       USE KEY722:       USE KEY
822:       implicit NONE723:       implicit NONE
823:       integer I, N, J, AN, NATOMS724:       integer I, N, J, AN, NATOMS
824: 725: 
825: 726: 
826:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),727:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
827:      Q grad(3*NATOMS), energy,STT,SA,SB,SC728:      Q grad(3*NATOMS), energy,STT,SA,SB,SC
828: 729: 
829:       integer C2,DINDEX,ii,jj,kk,k,l,iii,jjj730:       integer C1, C2, ii,jj,kk,k,l,iii,jjj
830:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 731:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
831:         INTEGER NTHREADS,MAXSEP732:         INTEGER NTHREADS
832: !$    INTEGER  OMP_GET_NUM_THREADS,733: !$    INTEGER  OMP_GET_NUM_THREADS,
833: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS734: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
834: 735: 
835:         LOGICAL NNCinc(NATOMS,MAXSEP)736:         INTEGER NNCinc(NATOMS,NATOMS)
836:         DOUBLE PRECISION dx,dy,dz737:         DOUBLE PRECISION dx,dy,dz
837:         integer tempN, alpha738:         integer tempN, alpha
838:         double precision Rdiff,Vfunc,Ffunc739:         double precision Rdiff,Vfunc,Ffunc
839:         double precision Rcut2,Rswitch2740:         double precision Rcut2,Rswitch2
840:         ! Ngrid is the number of atoms in that grid point741:         ! Ngrid is the number of atoms in that grid point
841:         ! grid is the array of atoms in each grid742:         ! grid is the array of atoms in each grid
842:         integer Ngrid,grid,maxgrid,maxpergrid743:         integer Ngrid,grid,maxgrid,maxpergrid
843:         ! number of atoms per grid, max744:         ! number of atoms per grid, max
844:         parameter (maxpergrid=25)745:         parameter (maxpergrid=50)
845:         ! dimensions of grid746:         ! dimensions of grid
846:         parameter (maxgrid=30)747:         parameter (maxgrid=15)
847:         dimension Ngrid(maxgrid,maxgrid,maxgrid),748:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
848:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)749:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
849:         integer MaxGridX,MaxGridY,MaxGridZ750:         integer MaxGridX,MaxGridY,MaxGridZ
850:         double precision gridsize,RD1751:         double precision gridsize,RD1
851:         double precision minX,minY,minZ,maxX,maxY,maxZ752:         double precision minX,minY,minZ,maxX,maxY,maxZ
852:         integer Xgrid,Ygrid,Zgrid,TID753:         integer Xgrid,Ygrid,Zgrid,TID
853:         Rdiff=NCcut-NCswitch754:         Rdiff=NCcut-NCswitch
854:         alpha=12755:         alpha=12
855: 756: 
 757: 
856:         GRIDSIZE=NCcut*1.01758:         GRIDSIZE=NCcut*1.01
857:         Rcut2=NCcut**2759:         Rcut2=NCcut**2
858:         Rswitch2=NCswitch**2760:         Rswitch2=NCswitch**2
859: 761: 
860:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 762:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
861:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )763:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
862: 764: 
863:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)765:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)
864: 766: 
865:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)767:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)
894:            endif796:            endif
895:         enddo797:         enddo
896: 798: 
897:         maxgridX=int((maxX-minX)/gridsize)+1799:         maxgridX=int((maxX-minX)/gridsize)+1
898:         maxgridY=int((maxY-minY)/gridsize)+1800:         maxgridY=int((maxY-minY)/gridsize)+1
899:         maxgridZ=int((maxZ-minZ)/gridsize)+1801:         maxgridZ=int((maxZ-minZ)/gridsize)+1
900: 802: 
901:         if(maxgridX .ge. maxgrid .or. 803:         if(maxgridX .ge. maxgrid .or. 
902:      Q  maxgridY .ge. maxgrid .or.804:      Q  maxgridY .ge. maxgrid .or.
903:      Q  maxgridZ .ge. maxgrid )then805:      Q  maxgridZ .ge. maxgrid )then
904:         write(*,*) 'ERROR: system got too big for grid searching...'806:         write(*,*) 'system got too big for grid searching...'
905: !        call abort807: !        call abort
906:         endif808:         endif
907: 809: 
908: 810: 
909: !!  Add a second grid that only includes atoms that are charged811: !!  Add a second grid that only includes atoms that are charged
910: 812: 
911:         do i=1,maxgrid813:         do i=1,maxgrid
912:          do j=1,maxgrid814:          do j=1,maxgrid
913:           do k=1,maxgrid815:           do k=1,maxgrid
914:                 Ngrid(i,j,k)=0816:                 Ngrid(i,j,k)=0
915:           enddo817:           enddo
916:          enddo818:          enddo
917:         enddo819:         enddo
918:         do i=1,NATOMS820:         do i=1,NATOMS
919: 821: 
920:                 Xgrid=int((X(i)-minX)/gridsize)+1822:                 Xgrid=int((X(i)-minX)/gridsize)+1
921:                 Ygrid=int((Y(i)-minY)/gridsize)+1823:                 Ygrid=int((Y(i)-minY)/gridsize)+1
922:                 Zgrid=int((Z(i)-minZ)/gridsize)+1824:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
923:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1825:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
924:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then826:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
925:                         write(*,*) 'ERROR: too many atoms in a grid'827:                         write(*,*) 'too many atoms in a grid'
926:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid828:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
927: !                        call abort829: !                        call abort
928:                 endif830:                 endif
929:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i831:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
930:         enddo832:         enddo
931: 833: 
932: 834: 
933:            tempN=0835:            tempN=0
934: 836: 
935: ! add a second loop that goes over charged atoms837: ! add a second loop that goes over charged atoms
936: 838: 
937: !$OMP PARALLEL839: !$OMP PARALLEL
938: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C2,DINDEX,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)  840: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)  
939: !$OMP& REDUCTION(+:ENERGY,grad)841: !$OMP& REDUCTION(+:ENERGY,grad)
940:         TID=0842:         TID=0
941:         NTHREADS=1;843:         NTHREADS=1;
942: !$      TID = OMP_GET_THREAD_NUM()844: !$      TID = OMP_GET_THREAD_NUM()
943: !$      NTHREADS = OMP_GET_NUM_THREADS()845: !$      NTHREADS = OMP_GET_NUM_THREADS()
944:         do i=1+TID,NATOMS,NTHREADS846:         do i=1+TID,NATOMS,NTHREADS
945: 847: 
946:           Xgrid=int((X(I)-minX)/gridsize)+1848: !!!        do i=1,NATOMS
947:           Ygrid=int((Y(I)-minY)/gridsize)+1849:            C1 = i
948:           Zgrid=int((Z(I)-minZ)/gridsize)+1850:                 Xgrid=int((X(i)-minX)/gridsize)+1
949:           do ii=XGRID-1,XGRID+1851:                 Ygrid=int((Y(i)-minY)/gridsize)+1
950:            do jj=YGRID-1,YGRID+1852:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
951:             do kk=ZGRID-1,ZGRID+1853:         
952:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.854:            do ii=XGRID-1,XGRID+1
 855:             do jj=YGRID-1,YGRID+1
 856:              do kk=ZGRID-1,ZGRID+1
 857:            if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
953:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then858:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
954:               do jjj=1,Ngrid(ii,jj,kk) 
955: 859: 
956:                C2=grid(ii,jj,kk,jjj)860:            do jjj=1,Ngrid(ii,jj,kk)
957:                DINDEX=C2-I861: 
958:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNCinc(C1,C12))) )then862:            C2=grid(ii,jj,kk,jjj)
959:                if(DINDEX .gt. MAXSEP .or. (DINDEX .gt. 0 .and. NNCinc(I,DINDEX)))then863:              if(NNCinc(i,C2) .eq. 1)then
960:                 dx = X(I) - X(C2)864: !!! SSS
961:                 dy = Y(I) - Y(C2)865:         dx = X(C1) - X(C2)
962:                 dz = Z(I) - Z(C2)866:         dy = Y(C1) - Y(C2)
963:                 r2 = dx**2 + dy**2 + dz**2867:         dz = Z(C1) - Z(C2)
964:                 if(r2 .le. Rcut2)then868:           r2 = dx**2 + dy**2 + dz**2
965:                  rm2 = 1/r2869:           if(r2 .le. Rcut2)then
966:                  rm14 = rm2**7870:              rm2 = 1/r2
967:                  energy=energy+STT*rm2**6+SC871:              rm14 = rm2**7
968:                  f_over_r=-STT*12.0*rm14872:                 energy=energy+STT*rm2**6+SC
969:                  RD1=sqrt(r2)-NCswitch873:                 f_over_r=-STT*12.0*rm14
970:                  if(r2 .gt. Rswitch2)then874:                 RD1=sqrt(r2)-NCswitch
971:                   f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)875:                 if(r2 .gt. Rswitch2)then
972:                   energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4876:                         f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
973:                  elseif(r2 .lt. Rswitch2)then877:                         energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4
 878:                 elseif(r2 .lt. Rswitch2)then
974:                 ! normal repulsive term879:                 ! normal repulsive term
975:                  else880:                 else
976:                 ! things should have fallen in one of the previous two...881:                 ! things should have fallen in one of the previous two...
977:                   write(*,*) 'something went wrong with switching function'882:                 write(*,*) 'something went wrong with switching function'
978: !                call abort883: !                call abort
979:                  endif884:                 endif
980:                  grad(I*3-2) = grad(I*3-2) + f_over_r * dx885:               grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx
981:                  grad(I*3-1) = grad(I*3-1) + f_over_r * dy886:               grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy
982:                  grad(I*3)   = grad(I*3)   + f_over_r * dz887:               grad(C1*3)   = grad(C1*3)   + f_over_r * dz
983:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx888:               grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
984:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy889:               grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
985:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz890:               grad(C2*3)   = grad(C2*3)   - f_over_r * dz
986:                  endif891: 
987:                endif892:             endif
988:               enddo893: 
989:              endif894: 
 895:            endif
 896:            enddo
 897:  
 898:           endif
 899: 
990:             enddo900:             enddo
991:            enddo901:            enddo
992:           enddo902:           enddo
993:          enddo903: 
 904:           enddo
994: !$OMP END PARALLEL905: !$OMP END PARALLEL
995: 906: 
996:       END907:       END
997: 908: 
998: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^909: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
999: 910: 
1000: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<911: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1001: ! FUNCTIONS for DH Calculations912: ! FUNCTIONS for DH Calculations
1002: !*****************************************************913: !*****************************************************
1003: 914: 
1014:       double precision kappa,r925:       double precision kappa,r
1015:       SBMDHVP=-kappa*exp(-kappa*r)/r-exp(-kappa*r)/(r**2)926:       SBMDHVP=-kappa*exp(-kappa*r)/r-exp(-kappa*r)/(r**2)
1016:       END927:       END
1017: !****************************************************928: !****************************************************
1018: 929: 
1019: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<930: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1020: !* SBMDHELEC computes the forces due to DH electrostatics        *931: !* SBMDHELEC computes the forces due to DH electrostatics        *
1021: !**********************************************************************932: !**********************************************************************
1022: 933: 
1023:       subroutine SBMDHELEC(x,y,z,grad, energy, natoms, 934:       subroutine SBMDHELEC(x,y,z,grad, energy, natoms, 
1024:      Q NNCinc,MAXSEP,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)935:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
1025: 936: 
1026:       USE KEY937:       USE KEY
1027:       implicit NONE938:       implicit NONE
1028:       integer I, N, J, AN, NATOMS,MAXSEP939:       integer I, N, J, AN, NATOMS
1029: 940: 
1030: 941: 
1031:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),942:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
1032:      Q grad(3*NATOMS), energy,STT,SA,SB,SC,943:      Q grad(3*NATOMS), energy,STT,SA,SB,SC,
1033:      Q SBMDHVP, SBMDHV944:      Q SBMDHVP, SBMDHV
1034:       integer C1, C2, C12,ii,jj,kk,k,l,iii,jjj945:       integer C1, C2, ii,jj,kk,k,l,iii,jjj
1035:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 946:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 
1036:       DOUBLE PRECISION A,B,C, D, COEF2, COEF3947:       DOUBLE PRECISION A,B,C, D, COEF2, COEF3
1037:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,948:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,
1038:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS949:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
1039: 950: 
1040:         LOGICAL NNCinc(NATOMS,MAXSEP)951:         INTEGER NNCinc(NATOMS,NATOMS)
1041:         DOUBLE PRECISION dx,dy,dz952:         DOUBLE PRECISION dx,dy,dz
1042:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)953:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)
1043:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS)954:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS)
1044:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut2955:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut2
1045:         ! Ngrid is the number of atoms in that grid point956:         ! Ngrid is the number of atoms in that grid point
1046:         ! grid is the array of atoms in each grid957:         ! grid is the array of atoms in each grid
1047:         integer Ngrid,grid,maxgrid,maxpergrid958:         integer Ngrid,grid,maxgrid,maxpergrid
1048:         ! number of atoms per grid, max959:         ! number of atoms per grid, max
1049:         parameter (maxpergrid=50)960:         parameter (maxpergrid=50)
1050:         ! dimensions of grid961:         ! dimensions of grid
1130: 1041: 
1131: 1042: 
1132:         do ii=2,SBMCHARGEON(1)1043:         do ii=2,SBMCHARGEON(1)
1133:                 i=SBMCHARGEON(ii)1044:                 i=SBMCHARGEON(ii)
1134: 1045: 
1135:                 Xgrid=int((X(i)-minX)/gridsize)+11046:                 Xgrid=int((X(i)-minX)/gridsize)+1
1136:                 Ygrid=int((Y(i)-minY)/gridsize)+11047:                 Ygrid=int((Y(i)-minY)/gridsize)+1
1137:                 Zgrid=int((Z(i)-minZ)/gridsize)+11048:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
1138:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+11049:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
1139:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then1050:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
1140:                         write(*,*) 'ERROR: too many atoms in a grid in DH'1051:                         write(*,*) 'too many atoms in a grid'
1141:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid1052:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
1142: !                        call abort1053: !                        call abort
1143:                 endif1054:                 endif
1144:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i1055:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1145:         enddo1056:         enddo
1146: 1057: 
1147: 1058: 
1148:            tempN=01059:            tempN=0
1149: 1060: 
1150: 1061: 
1151: !$OMP PARALLEL1062: !$OMP PARALLEL
1152: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,C12,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  1063: !$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)  
1153: !$OMP& REDUCTION(+:energy,grad)1064: !$OMP& REDUCTION(+:energy,grad)
1154:         TID=01065:         TID=0
1155:         NTHREADS=1;1066:         NTHREADS=1;
1156: !$      TID = OMP_GET_THREAD_NUM()1067: !$      TID = OMP_GET_THREAD_NUM()
1157: !$      NTHREADS = OMP_GET_NUM_THREADS()1068: !$      NTHREADS = OMP_GET_NUM_THREADS()
1158:         do iii=2+TID,SBMCHARGEON(1),NTHREADS1069:         do iii=2+TID,SBMCHARGEON(1),NTHREADS
1159: 1070: 
1160: 1071: 
1161:            C1=SBMCHARGEON(iii)1072:            C1=SBMCHARGEON(iii)
1162:                 Xgrid=int((X(C1)-minX)/gridsize)+11073:                 Xgrid=int((X(C1)-minX)/gridsize)+1
1165:         1076:         
1166:            do ii=XGRID-1,XGRID+11077:            do ii=XGRID-1,XGRID+1
1167:             do jj=YGRID-1,YGRID+11078:             do jj=YGRID-1,YGRID+1
1168:              do kk=ZGRID-1,ZGRID+11079:              do kk=ZGRID-1,ZGRID+1
1169:            if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1080:            if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1170:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1081:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1171: 1082: 
1172:            do jjj=1,Ngrid(ii,jj,kk)1083:            do jjj=1,Ngrid(ii,jj,kk)
1173: 1084: 
1174:            C2=grid(ii,jj,kk,jjj)1085:            C2=grid(ii,jj,kk,jjj)
 1086:            if(NNCinc(C1,C2) .eq. 1)then
1175: 1087: 
1176:            C12=C2-C11088:         dx = X(C1) - X(C2)
1177:            if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNCinc(C1,C12))) )then1089:         dy = Y(C1) - Y(C2)
1178: 1090:         dz = Z(C1) - Z(C2)
1179:  
1180:           dx = X(C1) - X(C2) 
1181:           dy = Y(C1) - Y(C2) 
1182:           dz = Z(C1) - Z(C2) 
1183:           r2 = dx**2 + dy**2 + dz**21091:           r2 = dx**2 + dy**2 + dz**2
1184:           if(r2 .le. DHcut2)then1092:           if(r2 .le. DHcut2)then
1185:            C2T=SBMCHARGE(C1)*SBMCHARGE(C2)1093:            C2T=SBMCHARGE(C1)*SBMCHARGE(C2)
1186:         ! add force evaluation now1094:         ! add force evaluation now
1187: 1095: 
1188:              rm2 = 1/r21096:              rm2 = 1/r2
1189:              r1=sqrt(r2)1097:              r1=sqrt(r2)
1190:                 energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)1098:                 energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)
1191: 1099: 
1192:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)1100:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0