hdiff output

r30675/SBM.f 2016-07-06 15:35:49.719416113 +0100 r30674/SBM.f 2016-07-06 15:35:50.067420804 +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=25000) 11:         parameter(NSBMMAX=25000)
 12:       INTEGER MAXSEP,MAXSEP1 12:       INTEGER MAXSEP,MAXSEP1
 13:       PARAMETER (MAXSEP1=30) 13:       PARAMETER (MAXSEP1=30)
 14:       LOGICAL NNCINC(NSBMMAX,MAXSEP1) 14:       LOGICAL NNCinc(NSBMMAX,MAXSEP1)
 15:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 15:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX),
 16:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF1(NSBMMAX*5), 16:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF1(NSBMMAX*5),
 17:      Q CONCOEF2(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 17:      Q CONCOEF2(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
 18:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2), 18:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2),
 19:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3), 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), 20:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5),
 21:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX) 21:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX)
 22:       INTEGER MAXEXCLUSIONS,NEXCLUDE 
 23:       PARAMETER (MAXEXCLUSIONS=20) 
 24:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS) 
 25:  
 26:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3) 22:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3)
 27:         integer CONTACTTYPE,SBMCHARGEON(NSBMMAX) 23:         integer CONTACTTYPE,SBMCHARGEON(NSBMMAX)
 28:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 24:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 29:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 25:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
 30:      Q SBMPRK,SBMPRX,CONCOEF1, CONCOEF2 26:      Q SBMPRK,SBMPRX,CONCOEF1, CONCOEF2
 31:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP, 27:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP,
 32:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC, 28:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC,
 33:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,NNEXL2,NEXCLUDE 29:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,MAXSEP
 34:         common /logical/ NNCINC 30:         COMMON /logical/ NNCinc
 35:         MAXSEP=MAXSEP1 31:         MAXSEP=MAXSEP1
 36:  
 37:         if(NATOMS.gt. NSBMMAX)then 32:         if(NATOMS.gt. NSBMMAX)then
 38:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 33:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 39:         STOP 34:         STOP
 40:         endif 35:         endif
 41:  36: 
 42:        if(.NOT.CALLED)then 37:        if(.NOT.CALLED)then
 43:         NEXCLUDE=MAXEXCLUSIONS 
 44:         write(*,*) 38:         write(*,*)
 45:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 39:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 46:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 40:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 47:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 41:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 48:         write(*,*) 42:         write(*,*)
 49:          43:         
 50:  44: 
 51:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP,  45:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP, 
 52:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2,  46:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2, 
 53:      Q NNEXL1,NNEXL2,NEXCLUDE,NNCINC,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT, 47:      Q NNCinc,MAXSEP,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,
 54:      Q NSBMMAX,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 48:      Q NSBMMAX,CONTACTTYPE,
  49:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 55:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or. 50:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or.
 56:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then 51:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then
 57:         write(*,*) 'increase array size' 52:         write(*,*) 'increase array size'
 58:         write(*,*) NBA, NTA,NPA,NC  53:         write(*,*) NBA, NTA,NPA,NC 
 59:         stop 54:         stop
 60:         endif  55:         endif 
 61:         CALLED=.TRUE. 56:         CALLED=.TRUE.
 62:         endIF 57:         endIF
 63:  
 64:  
 65: ! call the energy routine 58: ! call the energy routine
  59: 
 66:       call calc_energy_SBM(qo,natoms,GRAD,energy,Ib1,Ib2, 60:       call calc_energy_SBM(qo,natoms,GRAD,energy,Ib1,Ib2,
 67:      Q Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK, 61:      Q Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,
 68:      Q PHITYPE,PHISBM,IC,JC,CONCOEF1,CONCOEF2,  62:      Q PHITYPE,PHISBM,IC,JC,CONCOEF1,CONCOEF2, 
 69:      Q NNEXL1,NNEXL2,NEXCLUDE,NNCINC,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON, 63:      Q  NNCinc,MAXSEP,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,
 70:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 64:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
 71:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 65:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 72:       IF (STEST) THEN 66:       IF (STEST) THEN
 73:          PRINT '(A)','ERROR - second derivatives not available' 67:          PRINT '(A)','ERROR - second derivatives not available'
 74:          STOP 68:          STOP
 75:       ENDIF 69:       ENDIF
 76:       return 70:       return
 77:       end 71:       end
 78:  72: 
 79:  73: 
 80: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 74: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 81: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 75: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 82: !* startt then the velocities are assigned, otherwise, they are read   * 76: !* startt then the velocities are assigned, otherwise, they are read   *
 83: !* by selecting 2, or generated by selecting 3                         * 77: !* by selecting 2, or generated by selecting 3                         *
 84: !*********************************************************************** 78: !***********************************************************************
 85:  79: 
 86:       subroutine SBMINIT(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk, 80:       subroutine SBMinit(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,
 87:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2,  81:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2, 
 88:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNCINC,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch, 82:      Q NNCinc,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,
 89:      Q NCcut,STT,NSBMMAX,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 83:      Q NCcut,STT,NSBMMAX,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
 90:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 84:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 91:       USE KEY 85:       USE KEY
 92:       USE COMMONS, only: ATMASS 86:       USE COMMONS, only: ATMASS
 93:       implicit NONE 87:       implicit NONE
 94:  88: 
 95:         integer i,j,MaxCon,NATOMS,storage, dummy,  ANr,  89:         integer i,j,MaxCon,NNCmax,NATOMS,storage, dummy,  ANr, 
 96:      Q  Ib22, Ib21,IT1, JT1, KT1, IT2, JT2, KT2, IP1, JP1, 90:      Q  Ib22, Ib21,IT1, JT1, KT1, IT2, JT2, KT2, IP1, JP1,
 97:      Q KP1, LP1, IP2, JP2, KP2, 91:      Q KP1, LP1, IP2, JP2, KP2,
 98:      Q LP2, nBA1, nTA1, nPA1, nBA2, nTA2, nPA2,  ind1, ind2, ANt, 92:      Q LP2, nBA1, nTA1, nPA1, nBA2, nTA2, nPA2,  ind1, ind2, ANt,
 99:      Q  MDT1, MDT2, cl1,cl2,tempi,NSBMMAX,MAXEXCLUSIONS,NT,MAXSEP,MAXSEPTEMP 93:      Q  MDT1, MDT2, cl1, cl2,tempi,NSBMMAX,MAXSEP,MAXSEPTEMP
100:  94: 
101:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX),  95:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 
102:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), 96:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX),
103:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS), 97:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS),
104:      Q CONCOEF1(5*NSBMMAX),CONCOEF2(5*NSBMMAX)  98:      Q CONCOEF1(5*NSBMMAX),CONCOEF2(5*NSBMMAX) 
105:       INTEGER  Ib1(2*NSBMMAX), Ib2(2*NSBMMAX), IT(2*NSBMMAX), JT(2*NSBMMAX),  99:       INTEGER  Ib1(2*NSBMMAX), Ib2(2*NSBMMAX), IT(2*NSBMMAX), JT(2*NSBMMAX), 
106:      Q KT(2*NSBMMAX),IP(3*NSBMMAX), JP(3*NSBMMAX), KP(3*NSBMMAX),100:      Q KT(2*NSBMMAX),IP(3*NSBMMAX), JP(3*NSBMMAX), KP(3*NSBMMAX),
107:      Q LP(3*NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 101:      Q LP(3*NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 
108:      Q  NBA, NTA, NPA, NC,SBMCHARGEON(NATOMS)102:      Q  NBA, NTA, NPA, NC,SBMCHARGEON(NATOMS)
109:       INTEGER PHITYPE(3*NSBMMAX)103:       INTEGER PHITYPE(3*NSBMMAX)
110:       integer AA,BB,ANTEMP104:       integer AA,BB,ANTEMP
111:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS) 
112:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS) 
113:       DOUBLE PRECISION dx,dy,dz 
114:       double precision PI 
115:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut 
116:       integer CONTACTTYPE 
117:       character TMPCHAR 
118:       integer TMPINT,NUMCHARGES,nexc,I1,I2 
119:       double precision TMPREAL,concentration 
120:       LOGICAL NNCINC(NATOMS,MAXSEP)105:       LOGICAL NNCINC(NATOMS,MAXSEP)
 106:       DOUBLE PRECISION dx,dy,dz
 107:        double precision PI
 108:         DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
 109:         integer CONTACTTYPE
 110:         character TMPCHAR
 111:         integer TMPINT,NUMCHARGES,nexc,I1,I2
 112:         double precision TMPREAL,concentration
121: 113: 
122:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)114:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)
123:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)115:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
124:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS116: 
125:       DIMENSION EXCLUSIONS(NATOMS,MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS) 
126: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 117: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
127: !$OMP PARALLEL118: !$OMP PARALLEL
128: !$    NTHREADS = OMP_GET_NUM_THREADS()119: !$    NTHREADS = OMP_GET_NUM_THREADS()
129: !$      TID = OMP_GET_THREAD_NUM()120: !$      TID = OMP_GET_THREAD_NUM()
130: !$    if(TID .eq. 0)then 121: !$    if(TID .eq. 0)then 
131: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS122: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS
132: !$    endif123: !$    endif
133: !$OMP END PARALLEL124: !$OMP END PARALLEL
 125:       pi = 3.14159265358979323846264338327950288419716939937510
 126:       MAXSEPTEMP=0
 127:       NNCmax = NATOMS*NATOMS
 128:       MaxCon=NATOMS*5
134: 129: 
135: 130:       do i=1,NATOMS
136:       DO I=1,NATOMS131:         do j=1,30
137:         DO J=1,MAXSEP132:           NNCINC(i,j)=.TRUE.
138:           NNCINC(I,J)=.TRUE. 
139:         enddo133:         enddo
140:       enddo134:       enddo
141: 135: 
142: 136: 
143: 137: 
144:       pi = 3.14159265358979323846264338327950288419716939937510 
145:       MaxCon=NATOMS*5 
146:       MAXSEPTEMP=0 
147:       do i=1,NATOMS 
148:         NUMOFEXCLUSIONS(I)=0 
149:       enddo 
150:  
151: ! These lines read in the parameters.138: ! These lines read in the parameters.
152:         open(30, file='SBM.INP', status='old', access='sequential')139:         open(30, file='SBM.INP', status='old', access='sequential')
153:         write(*,*) 140:         write(*,*) 
154:         write(*,*) 'Reading the forcefield from SBM.INP'141:         write(*,*) 'Reading the forcefield from SBM.INP'
155:         write(*,*) 142:         write(*,*) 
156:         read(30,*)143:         read(30,*)
157:         read(30,*)144:         read(30,*)
158:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut145:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut
159:         read(30,*)146:         read(30,*)
160:         read(30,*) RSig, Reps, NCswitch, NCcut147:         read(30,*) RSig, Reps, NCswitch, NCcut
161:         write(*,*) 'RSig, Reps,NCswitch,NCcut'148:         write(*,*) 'RSig, Reps,NCswitch,NCcut'
162:         write(*,'(4F10.5)') RSig, Reps, NCswitch, NCcut149:         write(*,'(4F10.5)') RSig, Reps, NCswitch, NCcut
163:         write(*,*) 'Reading electrostatic parameters'150:         write(*,*) 'Reading electrostatic parameters'
164:         write(*,'(5F10.5)') PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut151:         write(*,'(5F10.5)') PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut
165:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is in units152:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is is units
166:         ! A^-1153:         ! A^-1
167:         KAPPA=0.32*sqrt(CONCENTRATION)154:         KAPPA=0.32*sqrt(CONCENTRATION)
168:         PREFACTOR=PREFACTOR/DC155:         PREFACTOR=PREFACTOR/DC
169:         read(30,*) ANtemp156:         read(30,*) ANtemp
170:         write(*,*) ANtemp, 'atoms'157:         write(*,*) ANtemp, 'atoms'
171:         NUMCHARGES=1158:         NUMCHARGES=1
172:         do i=1, ANtemp159:         do i=1, ANtemp
173:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)160:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)
174:           if(SBMCHARGE(i) .ne. 0)then161:           if(SBMCHARGE(i) .ne. 0)then
175:              NUMCHARGES=NUMCHARGES+1162:              NUMCHARGES=NUMCHARGES+1
182:         read(30,*) NC169:         read(30,*) NC
183:         write(*,*) NC, 'contacts'        170:         write(*,*) NC, 'contacts'        
184:         read(30,*) CONTACTTYPE171:         read(30,*) CONTACTTYPE
185:           if(NC .gt. MaxCon)then172:           if(NC .gt. MaxCon)then
186:              write(*,*) 'too many contacts'173:              write(*,*) 'too many contacts'
187:              STOP174:              STOP
188:           endif175:           endif
189: 176: 
190:         do i=1, NC177:         do i=1, NC
191:           read(30, *) IC(i), JC(i), Sigma, EpsC178:           read(30, *) IC(i), JC(i), Sigma, EpsC
192:           !! Add more contact types179: 
193:           CALL INCLUDEEXCLUSIONS(NATOMS,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP) 
194:           if(CONTACTTYPE .eq. 1)then180:           if(CONTACTTYPE .eq. 1)then
195:             ConCoef1(I)=EPSC*SIGMA**12-REPS*RSIG**12181:             ConCoef1(I)=EPSC*SIGMA**12-REPS*RSIG**12
196:             ConCoef2(I)=2*EPSC*SIGMA**6182:             ConCoef2(I)=2*EPSC*SIGMA**6
197:           elseif(CONTACTTYPE .eq. 2)then183:           elseif(CONTACTTYPE .eq. 2)then
198:             ConCoef1(I)=5*EPSC*SIGMA**12-REPS*RSIG**12184:             ConCoef1(I)=5*EPSC*SIGMA**12-REPS*RSIG**12
199:             ConCoef2(I)=6*EPSC*SIGMA**10185:             ConCoef2(I)=6*EPSC*SIGMA**10
200:           else186:           else
201:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.'187:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.'
202:           endif188:           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
203:         end do198:         end do
204: ! make a function to check and add 
205: 199: 
206:           read(30,*)200:           read(30,*)
207:           read(30,*) nBA201:           read(30,*) nBA
208:           write(*,*) nBA, 'bonds'202:           write(*,*) nBA, 'bonds'
209:         do i=1, nBA203:         do i=1, nBA
210:           read(30,*) Ib1(i), Ib2(i),Rb(i), bK(i)204:           read(30,*) Ib1(i), Ib2(i),Rb(i), bK(i)
211:           !!! ADD BONDS OF TYPE 6205:           AA=min(IB1(I),IB2(I))
212:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP)THEN206:           BB=max(IB1(I),IB2(I))
213:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS LESS THAN MAXSEP'207:           BB=BB-AA
214:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'208:           if(BB .gt. MAXSEPTEMP)then
215:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'209:             MAXSEPTEMP=BB
216:           ENDIF210:           endif
217:           CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          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: 
218:         end do216:         end do
219: 217: 
220:           read(30,*)218:           read(30,*)
221:           read(30,*) nTA219:           read(30,*) nTA
222:           write(*,*) nTA, 'angles'220:           write(*,*) nTA, 'angles'
223:         do i=1, nTA221:         do i=1, nTA
224:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)222:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)
225:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          223: 
226:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          224:           AA=min(IT(I),JT(I))
227:           CALL INCLUDEEXCLUSIONS(NATOMS,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          225:           BB=max(IT(I),JT(I))
228:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.226:           BB=BB-AA
229:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.227:           if(BB .gt. MAXSEP)then
230:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN228:                 write(*,*) 'ERROR: exclusions between atoms too far apart',AA,BB
231:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS LESS THAN MAXSEP'229:           endif
232:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'230:           NNCINC(AA, BB)=.FALSE.
233:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'231:           AA=min(IT(I),KT(I))
234:           ENDIF232:           BB=max(IT(I),KT(I))
235:         ENDDO233:           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:         enddo
236: 253: 
237:           read(30,*) 254:           read(30,*) 
238:           read(30,*) nPA255:           read(30,*) nPA
239:           write(*,*) nPA, 'dihedrals'256:           write(*,*) nPA, 'dihedrals'
240: 257: 
241: ! this reads in the dihedral angles and calculates the cosines and sines258: ! this reads in the dihedral angles and calculates the cosines and sines
242: ! in order to make the force and energy calculations easier, later.259: ! in order to make the force and energy calculations easier, later.
243:         do i=1, npA260:         do i=1, npA
244:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)261:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)
245:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),JP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          262:           AA=min(IP(I),JP(I))
246:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          263:           BB=max(IP(I),JP(I))
247:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          264:           BB=BB-AA
248:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          265:           if(BB .gt. MAXSEPTEMP)then
249:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          266:            MAXSEPTEMP=BB
250:           CALL INCLUDEEXCLUSIONS(NATOMS,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          267:           endif
251:  
252:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR. 
253:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR. 
254:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR. 
255:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR. 
256:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR. 
257:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN 
258:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS LESS THAN MAXSEP' 
259:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.' 
260:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.' 
261:           ENDIF 
262:         ENDDO  
263: 268: 
264:         read(30,*)269:           if(BB .gt. MAXSEP)then
265:         read(30,*) nexc270:                 write(*,*) 'ERROR: exclusions between atoms too far apart.',AA,BB
266:         write(*,*) nexc, 'exculusions'271:           endif
 272:           NNCINC(AA, BB)=.FALSE.
 273:           AA=min(IP(I),KP(I))
 274:           BB=max(IP(I),KP(I))
 275:           BB=BB-AA
 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 DO
 329: 
 330:           read(30,*)
 331:           read(30,*) nexc
 332:           write(*,*) nexc, 'exculusions'
 333: 
 334:         do i=1, nexc
 335:           read(30,*) I1, I2
 336:           AA=min(I1,I2)
 337:           BB=max(I1,I2)
 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 do
267: 348: 
268:       do i=1, nexc 
269:         read(30,*) I1, I2 
270:         CALL INCLUDEEXCLUSIONS(NATOMS,I1,I2,EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)           
271:       enddo 
272: 349: 
273: ! now that we know what to exclude, let's make a static list of pairs to 
274: ! exclude 
275: ! At this point, we will change the meaning of maxexclusions, which is 
276: ! actually NEXCLUDE 
277:         NT=MAXEXCLUSIONS*NATOMS       
278:         MAXEXCLUSIONS=0   
279:         DO I=1,NATOMS 
280:           DO J=1,NUMOFEXCLUSIONS(I) 
281:             MAXEXCLUSIONS=MAXEXCLUSIONS+1 
282:             IF(MAXEXCLUSIONS .GT. NT)THEN 
283:               write(*,*) 'ERROR: TOO MANY EXCLUSIONS, IN TOTAL.' 
284:             ENDIF 
285:             NNEXL1(MAXEXCLUSIONS)=I 
286:             NNEXL2(MAXEXCLUSIONS)=EXCLUSIONS(I,J) 
287:           ENDDO 
288:         ENDDO 
289:         STT=reps*rsig**12350:         STT=reps*rsig**12
290: 351: 
291:         do I=1,NATOMS352:         do I=1,NATOMS
292:           TARR(I)=0353:           TARR(I)=0
293:         enddo354:         enddo
294: 355: 
295: ! read in position restraints356: ! read in position restraints
296:        read(30,*) 357:        read(30,*) 
297:        read(30,*) SBMPRN358:        read(30,*) SBMPRN
298:        write(*,*) SBMPRN, 'position restraints'359:        write(*,*) SBMPRN, 'position restraints'
304:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)365:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
305: !                call abort366: !                call abort
306:             endif367:             endif
307:             TARR(SBMPRI(I))=1368:             TARR(SBMPRI(I))=1
308:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  369:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  
309:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'370:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'
310:         endif371:         endif
311:        enddo 372:        enddo 
312:         MAXSEP=MAXSEPTEMP373:         MAXSEP=MAXSEPTEMP
313:        close(30)374:        close(30)
314:                 I=1381 
315:                 write(*,*) I,CONCOEF1(I),CONCOEF2(I) 
316:        end375:        end
317: 376: 
318: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^377: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
319: 378: 
320: 379: 
321: !************************************************************************** 
322: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in 
323: ! the model 
324: !************************************************************************** 
325:  
326:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE, 
327:      Q NNCINC,MAXSEP,MAXSEPTEMP) 
328:       IMPLICIT NONE 
329:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,MAXSEPTEMP 
330:       DIMENSION EXCLIST(NATOMS,MAXEXCLUDE),NEXCLUDE(NATOMS) 
331:       LOGICAL NNCINC(NATOMS,MAXSEP) 
332:       INTEGER I,N,M,AA,BB,AB 
333:       LOGICAL INCLUDED 
334:  
335: ! If the atoms are within MAXSEP of each other, then keep track via 
336: ! NNCINC 
337:       AA=min(ATOM1,ATOM2) 
338:       BB=max(ATOM1,ATOM2) 
339:       AB=BB-AA 
340:       IF(AB .le. MAXSEP)THEN 
341:         NNCINC(AA, AB)=.FALSE. 
342:         IF(AB .gt. MAXSEPTEMP)THEN 
343:           MAXSEPTEMP=AB 
344:         ENDIF 
345:       ELSE 
346:         INCLUDED = .FALSE. 
347:         N=MAX(ATOM1,ATOM2) 
348:         M=MIN(ATOM1,ATOM2) 
349:         DO I =1,NEXCLUDE(M) 
350:           IF(EXCLIST(M,I) .eq. N)THEN 
351:             INCLUDED = .TRUE. 
352:           ENDIF 
353:         ENDDO  
354:    
355:         if(.NOT. INCLUDED)THEN 
356:           if(NEXCLUDE(M) .EQ. MAXEXCLUDE)THEN 
357:             write(*,*) 'ERROR: TOO MANY EXCLUSIONS WITH ATOM',M 
358:           ENDIF 
359:           NEXCLUDE(M)=NEXCLUDE(M)+1 
360:           EXCLIST(M,NEXCLUDE(M))=N 
361:         ENDIF 
362:       ENDIF 
363:       END 
364: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
365:  
366:  
367: C380: C
368: C Calculate the Forces and energies381: C Calculate the Forces and energies
369: C382: C
370:       subroutine CALC_ENERGY_SBM(qo,natoms,GRAD, energy,Ib1, Ib2,383:       subroutine calc_energy_SBM(qo,natoms,GRAD, energy,Ib1, Ib2,
371:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,384:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,
372:      Q PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2, 385:      Q PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2, 
373:      Q NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE, SBMCHARGEON,NCswitch,386:      Q NNCinc,MAXSEP,NBA, NTA, NPA, NC,SBMCHARGE, SBMCHARGEON,NCswitch,
374:      Q STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,387:      Q STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
375:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)388:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
376: 389: 
377:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP390:       INTEGER I, J, NATOMS,NBA, NTA, NPA, NC,MAXSEP
378: 391: 
379:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY392:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
380:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)393:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
381:       LOGICAL NNCINC(NATOMS,MAXSEP) 
382: 394: 
383:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 395:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
384:      Q PHISBM(NPA),CONCOEF1(NC),CONCOEF2(NC),NCswitch,NCcut,STT, 396:      Q PHISBM(NPA),CONCOEF1(NC),CONCOEF2(NC),NCswitch,NCcut,STT, 
385:      Q SBMCHARGE(NATOMS)397:      Q SBMCHARGE(NATOMS)
386:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 398:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 
387:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),SBMCHARGEON(NATOMS),399:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),SBMCHARGEON(NATOMS),
388:      Q PHITYPE(NPA),CONTACTTYPE400:      Q PHITYPE(NPA),CONTACTTYPE
389:         INTEGER NNEXL1(NEXCLUSIONS)401:         LOGICAL NNCinc(NATOMS,MAXSEP)
390:         INTEGER NNEXL2(NEXCLUSIONS) 
391:       DOUBLE PRECISION dx,dy,dz402:       DOUBLE PRECISION dx,dy,dz
392:       INTEGER SBMPRN, SBMPRI(NATOMS)403:       INTEGER SBMPRN, SBMPRI(NATOMS)
393:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)404:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
394:       do i = 1, natoms405:       do i = 1, natoms
395:         j = (i-1)*3406:         j = (i-1)*3
396:         x(i) = qo(j+1)407:         x(i) = qo(j+1)
397:         y(i) = qo(j+2)408:         y(i) = qo(j+2)
398:         z(i) = qo(j+3)409:         z(i) = qo(j+3)
399:         grad(j+1) = 0.0410:         grad(j+1) = 0.0
400:         grad(j+2) = 0.0411:         grad(j+2) = 0.0
402:       enddo413:       enddo
403: 414: 
404:       energy = 0.0415:       energy = 0.0
405:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)416:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
406:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)417:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
407:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,418:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
408:      Q PHITYPE,PHISBM,NPA)419:      Q PHITYPE,PHISBM,NPA)
409:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC, 420:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC, 
410:      Q CONCOEF1,CONCOEF2,NC,CONTACTTYPE)421:      Q CONCOEF1,CONCOEF2,NC,CONTACTTYPE)
411:       call SBMNonContacts(x,y,z,grad, energy, natoms, 422:       call SBMNonContacts(x,y,z,grad, energy, natoms, 
412:      Q NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,NCswitch,NCcut,STT)423:      Q NNCinc,MAXSEP,NCswitch,NCcut,STT)
413: !      call SBMDHELEC(x,y,z,grad, energy, natoms, 424:       call SBMDHELEC(x,y,z,grad, energy, natoms, 
414: !     Q NNEXL1,NNEXL2,NEXCLUSIONS,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)425:      Q NNCinc,MAXSEP,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
415:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)426:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
416: 427: 
417:       end428:       end
418: 429: 
419: 430: 
420: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<431: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
421: !* SBMBonds  computes the hookean force and energy between chosen atoms *432: !* SBMBonds  computes the hookean force and energy between chosen atoms *
422: !***********************************************************************433: !***********************************************************************
423: 434: 
424:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)435:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)
731: 742: 
732: ! type 1 is 6-12743: ! type 1 is 6-12
733: ! type 2 is 12-12744: ! type 2 is 12-12
734: ! implement type 3 as gaussian745: ! implement type 3 as gaussian
735: ! implement type 4 as dual-basin gaussian746: ! implement type 4 as dual-basin gaussian
736: ! question: Should type be part of the interaction, or should they be organized747: ! question: Should type be part of the interaction, or should they be organized
737: ! into separate lists?748: ! into separate lists?
738: 749: 
739: ! type 1 is 6-12 interaction750: ! type 1 is 6-12 interaction
740:         if(CONTACTTYPE .eq. 1)then751:         if(CONTACTTYPE .eq. 1)then
741: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R) REDUCTION(+:ENERGY,grad)752: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R)REDUCTION(+:ENERGY,grad)
742: !$OMP DO753: !$OMP DO
743:           do i=1, NC754:           do i=1, NC
744:             C1 = IC(i)755:             C1 = IC(i)
745:             C2 = JC(i)756:             C2 = JC(i)
746:             dx = X(C1) - X(C2)757:             dx = X(C1) - X(C2)
747:             dy = Y(C1) - Y(C2)758:             dy = Y(C1) - Y(C2)
748:             dz = Z(C1) - Z(C2)759:             dz = Z(C1) - Z(C2)
749:             r2 = dx**2 + dy**2 + dz**2760:             r2 = dx**2 + dy**2 + dz**2
750:             rm2 = 1.0/r2761:             rm2 = 1.0/r2
751:           !  rm2 = rm2*sigma(i)762:           !  rm2 = rm2*sigma(i)
752:             rm6 = rm2**3763:             rm6 = rm2**3
753:             !ENERGY = ENERGY + epsC(i)*rm6*(rm6-2.0)764:             !ENERGY = ENERGY + epsC(i)*rm6*(rm6-2.0)
754:             ENERGY = ENERGY + rm6*(CONCOEF1(I)*rm6-CONCOEF2(I))765:             ENERGY = ENERGY + rm6*(CONCOEF1(I)*rm6-CONCOEF2(I))
755:             f_over_r = -12.0*rm6*(CONCOEF1(I)*rm6-CONCOEF2(I)/2.0)*rm2766:             f_over_r = -12.0*rm6*(CONCOEF1(I)*rm6-CONCOEF2(I)/2.0)*rm2
756:  
757:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx767:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
758:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy768:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
759:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz769:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz
760:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx770:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
761:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy771:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
762:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz772:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz
763:           enddo773:           enddo
764: !$OMP END DO774: !$OMP END DO
765: !$OMP END PARALLEL775: !$OMP END PARALLEL
766: ! type 2 is 10-12 interaction776: ! type 2 is 10-12 interaction
767:         elseif(CONTACTTYPE .eq. 2)then777:         elseif(CONTACTTYPE .eq. 2)then
768: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R) REDUCTION(+:ENERGY,grad)778: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R)REDUCTION(+:ENERGY,grad)
769: !$OMP DO779: !$OMP DO
770:           do i=1, NC780:           do i=1, NC
771:              C1 = IC(i)781:              C1 = IC(i)
772:              C2 = JC(i)782:              C2 = JC(i)
773:              dx = X(C1) - X(C2)783:              dx = X(C1) - X(C2)
774:              dy = Y(C1) - Y(C2)784:              dy = Y(C1) - Y(C2)
775:              dz = Z(C1) - Z(C2)785:              dz = Z(C1) - Z(C2)
776:              r2 = dx**2 + dy**2 + dz**2786:              r2 = dx**2 + dy**2 + dz**2
777:              rm2 = 1.0/r2787:              rm2 = 1.0/r2
778: !             rm2 = rm2*sigma(i)788: !             rm2 = rm2*sigma(i)
798: 808: 
799:       end809:       end
800: 810: 
801: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^811: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
802: 812: 
803: 813: 
804: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<814: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
805: !* SBMNonContacts computes the forces due to non native contacts       *815: !* SBMNonContacts computes the forces due to non native contacts       *
806: !**********************************************************************816: !**********************************************************************
807: 817: 
808:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,NNEXL1,NNEXL2,818:       subroutine SBMnoncontacts(x,y,z,grad, energy, 
809:      Q  NEXCLUSIONS,NNCINC,MAXSEP,NCswitch,NCcut,STT)819:      Q NATOMS,NNCinc,MAXSEP,NCswitch,NCcut,STT)
810:         USE KEY820:       USE KEY
811:         implicit NONE821:       implicit NONE
812:         integer I, N, J, AN, NATOMS822:       integer I, N, J, AN, NATOMS
 823: 
813: 824: 
 825:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
 826:      Q grad(3*NATOMS), energy,STT,SA,SB,SC
814: 827: 
815:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),828:       integer C2,DINDEX,ii,jj,kk,k,l,iii,jjj
816:      Q  grad(3*NATOMS), energy,STT,SA,SB,SC829:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
817: 830:         INTEGER NTHREADS,MAXSEP
818:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj 
819:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut  
820:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP 
821:         LOGICAL NNCINC(NATOMS,MAXSEP) 
822: !$    INTEGER  OMP_GET_NUM_THREADS,831: !$    INTEGER  OMP_GET_NUM_THREADS,
823: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS832: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
824: 833: 
825:         INTEGER NNEXL1(NEXCLUSIONS)834:         LOGICAL NNCinc(NATOMS,MAXSEP)
826:         INTEGER NNEXL2(NEXCLUSIONS) 
827:         DOUBLE PRECISION dx,dy,dz835:         DOUBLE PRECISION dx,dy,dz
828:         integer tempN, alpha836:         integer tempN, alpha
829:         double precision Rdiff,Vfunc,Ffunc837:         double precision Rdiff,Vfunc,Ffunc
830:         double precision Rcut2,Rswitch2838:         double precision Rcut2,Rswitch2
831:         ! Ngrid is the number of atoms in that grid point839:         ! Ngrid is the number of atoms in that grid point
832:         ! grid is the array of atoms in each grid840:         ! grid is the array of atoms in each grid
833:         integer Ngrid,grid,maxgrid,maxpergrid841:         integer Ngrid,grid,maxgrid,maxpergrid
834:         ! number of atoms per grid, max842:         ! number of atoms per grid, max
835:         parameter (maxpergrid=50)843:         parameter (maxpergrid=50)
836:         ! dimensions of grid844:         ! dimensions of grid
919:                 endif927:                 endif
920:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i928:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
921:         enddo929:         enddo
922: 930: 
923: 931: 
924:            tempN=0932:            tempN=0
925: 933: 
926: ! add a second loop that goes over charged atoms934: ! add a second loop that goes over charged atoms
927: 935: 
928: !$OMP PARALLEL936: !$OMP PARALLEL
929: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)  937: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C2,DINDEX,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)  
930: !$OMP& REDUCTION(+:ENERGY,grad)938: !$OMP& REDUCTION(+:ENERGY,grad)
931:         TID=0939:         TID=0
932:         NTHREADS=1;940:         NTHREADS=1;
933: !$      TID = OMP_GET_THREAD_NUM()941: !$      TID = OMP_GET_THREAD_NUM()
934: !$      NTHREADS = OMP_GET_NUM_THREADS()942: !$      NTHREADS = OMP_GET_NUM_THREADS()
935:         do i=1+TID,NATOMS,NTHREADS943:         do i=1+TID,NATOMS,NTHREADS
936: 944: 
937:           Xgrid=int((X(I)-minX)/gridsize)+1945:           Xgrid=int((X(I)-minX)/gridsize)+1
938:           Ygrid=int((Y(I)-minY)/gridsize)+1946:           Ygrid=int((Y(I)-minY)/gridsize)+1
939:           Zgrid=int((Z(I)-minZ)/gridsize)+1947:           Zgrid=int((Z(I)-minZ)/gridsize)+1
940:           do ii=XGRID-1,XGRID+1948:           do ii=XGRID-1,XGRID+1
941:            do jj=YGRID-1,YGRID+1949:            do jj=YGRID-1,YGRID+1
942:             do kk=ZGRID-1,ZGRID+1950:             do kk=ZGRID-1,ZGRID+1
943:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.951:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
944:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then952:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
945:               do jjj=1,Ngrid(ii,jj,kk)953:               do jjj=1,Ngrid(ii,jj,kk)
946: 954: 
947:                C2=grid(ii,jj,kk,jjj)955:                C2=grid(ii,jj,kk,jjj)
948:                DINDEX=C2-I956:                DINDEX=C2-I
949:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNEXL(C1,C12))) )then957:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNCinc(C1,C12))) )then
950:            !if(DINDEX .gt. MAXSEP .or. (DINDEX .gt. 0 .and. NNEXL(I,DINDEX)))then958:                if(DINDEX .gt. MAXSEP .or. (DINDEX .gt. 0 .and. NNCinc(I,DINDEX)))then
951:                if(DINDEX .gt. 0 .AND. ( DINDEX .GT. MAXSEP .OR. (DINDEX .LE. MAXSEP .AND. NNCINC(I,DINDEX) ) ) )then 
952:                 dx = X(I) - X(C2)959:                 dx = X(I) - X(C2)
953:                 dy = Y(I) - Y(C2)960:                 dy = Y(I) - Y(C2)
954:                 dz = Z(I) - Z(C2)961:                 dz = Z(I) - Z(C2)
955:                 r2 = dx**2 + dy**2 + dz**2962:                 r2 = dx**2 + dy**2 + dz**2
956:                 if(r2 .le. Rcut2)then963:                 if(r2 .le. Rcut2)then
957:                  rm2 = 1/r2964:                  rm2 = 1/r2
958:                  rm14 = rm2**7965:                  rm14 = rm2**7
959:                  energy=energy+STT*rm2**6+SC966:                  energy=energy+STT*rm2**6+SC
960:                  f_over_r=-STT*12.0*rm14967:                  f_over_r=-STT*12.0*rm14
961:                  RD1=sqrt(r2)-NCswitch968:                  RD1=sqrt(r2)-NCswitch
962:                  if(r2 .gt. Rswitch2)then969:                  if(r2 .gt. Rswitch2)then
963:                   f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)970:                   f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
964:                   energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4971:                   energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4
965:                  elseif(r2 .lt. Rswitch2)then972:                  elseif(r2 .lt. Rswitch2)then
966:                 ! normal repulsive term973:                 ! normal repulsive term
967:                  else974:                  else
968:                 ! things should have fallen in one of the previous two...975:                 ! things should have fallen in one of the previous two...
969:                   write(*,*) 'something went wrong with switching function'976:                   write(*,*) 'something went wrong with switching function'
 977: !                call abort
970:                  endif978:                  endif
971:                  grad(I*3-2) = grad(I*3-2) + f_over_r * dx979:                  grad(I*3-2) = grad(I*3-2) + f_over_r * dx
972:                  grad(I*3-1) = grad(I*3-1) + f_over_r * dy980:                  grad(I*3-1) = grad(I*3-1) + f_over_r * dy
973:                  grad(I*3)   = grad(I*3)   + f_over_r * dz981:                  grad(I*3)   = grad(I*3)   + f_over_r * dz
974:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx982:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
975:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy983:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
976:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz984:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz
977:                  endif985:                  endif
978:                endif986:                endif
979:               enddo987:               enddo
980:              endif988:              endif
981:             enddo989:             enddo
982:            enddo990:            enddo
983:           enddo991:           enddo
984:          enddo992:          enddo
985: !$OMP END PARALLEL993: !$OMP END PARALLEL
986: 994: 
987:        DO I=1,NEXCLUSIONS 
988:          C1=NNEXL1(I)  
989:          C2=NNEXL2(I) 
990:          dx = X(C1) - X(C2) 
991:          dy = Y(C1) - Y(C2) 
992:          dz = Z(C1) - Z(C2) 
993:          r2 = dx**2 + dy**2 + dz**2 
994:          if(r2 .le. Rcut2)then 
995:           rm2 = 1/r2 
996:           rm14 = rm2**7 
997:           energy=energy-STT*rm2**6-SC 
998:           f_over_r=-STT*12.0*rm14 
999:           RD1=sqrt(r2)-NCswitch 
1000:           if(r2 .gt. Rswitch2)then 
1001:            f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2) 
1002:            energy=energy-SA/3.0*RD1**3-SB/4.0*RD1**4 
1003:           elseif(r2 .lt. Rswitch2)then 
1004:          ! normal repulsive term 
1005:           else 
1006:          ! things should have fallen in one of the previous two... 
1007:            write(*,*) 'something went wrong with switching function' 
1008:           endif 
1009:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx 
1010:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy 
1011:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz 
1012:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx 
1013:           grad(C2*3-1) = grad(C2*3-1) + f_over_r * dy 
1014:           grad(C2*3)   = grad(C2*3)   + f_over_r * dz 
1015:          endif 
1016:       ENDDO 
1017:  
1018:       END995:       END
1019: 996: 
1020: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^997: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
1021: 998: 
1022: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<999: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1023: ! FUNCTIONS for DH Calculations1000: ! FUNCTIONS for DH Calculations
1024: !*****************************************************1001: !*****************************************************
1025: 1002: 
1026:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)1003:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)
1027:       USE KEY1004:       USE KEY
1036:       double precision kappa,r1013:       double precision kappa,r
1037:       SBMDHVP=-kappa*exp(-kappa*r)/r-exp(-kappa*r)/(r**2)1014:       SBMDHVP=-kappa*exp(-kappa*r)/r-exp(-kappa*r)/(r**2)
1038:       END1015:       END
1039: !****************************************************1016: !****************************************************
1040: 1017: 
1041: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1018: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1042: !* SBMDHELEC computes the forces due to DH electrostatics        *1019: !* SBMDHELEC computes the forces due to DH electrostatics        *
1043: !**********************************************************************1020: !**********************************************************************
1044: 1021: 
1045:       subroutine SBMDHELEC(x,y,z,grad, energy, natoms, 1022:       subroutine SBMDHELEC(x,y,z,grad, energy, natoms, 
1046:      Q NNEXL1,NNEXL2,NEXCLUSIONS,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)1023:      Q NNCinc,MAXSEP,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
1047: 1024: 
1048:       USE KEY1025:       USE KEY
1049:       implicit NONE1026:       implicit NONE
1050:       integer I, N, J, AN, NATOMS1027:       integer I, N, J, AN, NATOMS,MAXSEP
1051: 1028: 
1052: 1029: 
1053:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),1030:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
1054:      Q grad(3*NATOMS), energy,STT,SA,SB,SC,1031:      Q grad(3*NATOMS), energy,STT,SA,SB,SC,
1055:      Q SBMDHVP, SBMDHV1032:      Q SBMDHVP, SBMDHV
1056:       integer C1, C2, C12,ii,jj,kk,k,l,iii,jjj1033:       integer C1, C2, C12,ii,jj,kk,k,l,iii,jjj
1057:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 1034:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 
1058:       DOUBLE PRECISION A,B,C, D, COEF2, COEF31035:       DOUBLE PRECISION A,B,C, D, COEF2, COEF3
1059:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,1036:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,
1060:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS1037:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
1061:         INTEGER NEXCLUSIONS1038: 
1062:         INTEGER NNEXL1(NEXCLUSIONS)1039:         LOGICAL NNCinc(NATOMS,MAXSEP)
1063:         INTEGER NNEXL2(NEXCLUSIONS) 
1064:         DOUBLE PRECISION dx,dy,dz1040:         DOUBLE PRECISION dx,dy,dz
1065:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)1041:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)
1066:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS)1042:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS)
1067:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut21043:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut2
1068:         ! Ngrid is the number of atoms in that grid point1044:         ! Ngrid is the number of atoms in that grid point
1069:         ! grid is the array of atoms in each grid1045:         ! grid is the array of atoms in each grid
1070:         integer Ngrid,grid,maxgrid,maxpergrid1046:         integer Ngrid,grid,maxgrid,maxpergrid
1071:         ! number of atoms per grid, max1047:         ! number of atoms per grid, max
1072:         parameter (maxpergrid=50)1048:         parameter (maxpergrid=50)
1073:         ! dimensions of grid1049:         ! dimensions of grid
1171:            tempN=01147:            tempN=0
1172: 1148: 
1173: 1149: 
1174: !$OMP PARALLEL1150: !$OMP PARALLEL
1175: !$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)  1151: !$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)  
1176: !$OMP& REDUCTION(+:energy,grad)1152: !$OMP& REDUCTION(+:energy,grad)
1177:         TID=01153:         TID=0
1178:         NTHREADS=1;1154:         NTHREADS=1;
1179: !$      TID = OMP_GET_THREAD_NUM()1155: !$      TID = OMP_GET_THREAD_NUM()
1180: !$      NTHREADS = OMP_GET_NUM_THREADS()1156: !$      NTHREADS = OMP_GET_NUM_THREADS()
1181:          do iii=2+TID,SBMCHARGEON(1),NTHREADS1157:         do iii=2+TID,SBMCHARGEON(1),NTHREADS
1182:           C1=SBMCHARGEON(iii)1158: 
1183:           Xgrid=int((X(C1)-minX)/gridsize)+11159: 
1184:           Ygrid=int((Y(C1)-minY)/gridsize)+11160:            C1=SBMCHARGEON(iii)
1185:           Zgrid=int((Z(C1)-minZ)/gridsize)+11161:                 Xgrid=int((X(C1)-minX)/gridsize)+1
 1162:                 Ygrid=int((Y(C1)-minY)/gridsize)+1
 1163:                 Zgrid=int((Z(C1)-minZ)/gridsize)+1
1186:         1164:         
1187:           do ii=XGRID-1,XGRID+11165:            do ii=XGRID-1,XGRID+1
1188:            do jj=YGRID-1,YGRID+11166:             do jj=YGRID-1,YGRID+1
1189:             do kk=ZGRID-1,ZGRID+11167:              do kk=ZGRID-1,ZGRID+1
1190:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1168:            if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1191:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1169:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1192:               do jjj=1,Ngrid(ii,jj,kk)1170: 
1193:                C2=grid(ii,jj,kk,jjj)1171:            do jjj=1,Ngrid(ii,jj,kk)
1194:            !C12=C2-C11172: 
1195:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNEXL(C1,C12))) )then1173:            C2=grid(ii,jj,kk,jjj)
1196:                if(C2-C1 .gt. 0)then1174: 
1197:                 dx = X(C1) - X(C2)1175:            C12=C2-C1
1198:                 dy = Y(C1) - Y(C2)1176:            if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNCinc(C1,C12))) )then
1199:                 dz = Z(C1) - Z(C2)1177: 
1200:                 r2 = dx**2 + dy**2 + dz**21178: 
1201:                 if(r2 .le. DHcut2)then1179:           dx = X(C1) - X(C2)
1202:                  C2T=SBMCHARGE(C1)*SBMCHARGE(C2)1180:           dy = Y(C1) - Y(C2)
 1181:           dz = Z(C1) - Z(C2)
 1182:           r2 = dx**2 + dy**2 + dz**2
 1183:           if(r2 .le. DHcut2)then
 1184:            C2T=SBMCHARGE(C1)*SBMCHARGE(C2)
1203:         ! add force evaluation now1185:         ! add force evaluation now
1204:                  rm2 = 1/r21186: 
1205:                  r1=sqrt(r2)1187:              rm2 = 1/r2
1206:                  energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)1188:              r1=sqrt(r2)
1207:                  f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)1189:                 energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)
1208:                  if(r2 .gt. DHswitch2)then1190: 
1209:                   RD1=r1-DHswitch1191:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)
1210:                   f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))1192:                 if(r2 .gt. DHswitch2)then
1211:                   energy=energy+COEF2*RD1**2+COEF3*RD1**31193:                         RD1=r1-DHswitch
1212:                  elseif(r2 .lt. DHswitch2)then1194:                         f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))
 1195:                         energy=energy+COEF2*RD1**2+COEF3*RD1**3
 1196:                 elseif(r2 .lt. DHswitch2)then
1213:                 ! normal DH term1197:                 ! normal DH term
1214:                  else1198:                 else
1215:                 ! things should have fallen in one of the previous two...1199:                 ! things should have fallen in one of the previous two...
1216:                   write(*,*) 'something went wrong with DH switching function'1200:                 write(*,*) 'something went wrong with DH switching function'
1217: !                call abort1201: !                call abort
1218:                  endif1202:                 endif
1219:                  f_over_r=f_over_r*sqrt(rm2)1203:                 f_over_r=f_over_r*sqrt(rm2)
1220:                  grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx1204:               grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx
1221:                  grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy1205:               grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy
1222:                  grad(C1*3)   = grad(C1*3)   + f_over_r * dz1206:               grad(C1*3)   = grad(C1*3)   + f_over_r * dz
1223:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx1207:               grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
1224:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy1208:               grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
1225:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz1209:               grad(C2*3)   = grad(C2*3)   - f_over_r * dz
 1210: 
 1211:             endif
 1212: 
 1213:            endif
 1214:            enddo
 1215:  
 1216:           endif
1226: 1217: 
1227:                   endif 
1228:                 endif 
1229:                enddo 
1230:                   endif 
1231:              enddo 
1232:             enddo1218:             enddo
1233:            enddo1219:            enddo
1234:           enddo1220:           enddo
1235: 1221: 
1236: !!! Add routine to subtract exclusion interactions1222:           enddo
1237: !$OMP DO 
1238:        DO I=1,NEXCLUSIONS 
1239:          C1=NNEXL1(I)  
1240:          C2=NNEXL2(I) 
1241:          dx = X(C1) - X(C2) 
1242:          dy = Y(C1) - Y(C2) 
1243:          dz = Z(C1) - Z(C2) 
1244:          r2 = dx**2 + dy**2 + dz**2 
1245:          if(r2 .le. DHcut2)then 
1246:           C2T=SBMCHARGE(C1)*SBMCHARGE(C2) 
1247:  ! add force evaluation now 
1248:           rm2 = 1/r2 
1249:           r1=sqrt(r2) 
1250:           energy=energy-PREFACTOR*C2T*SBMDHV(kappa,r1) 
1251:           f_over_r=-PREFACTOR*C2T*SBMDHVP(kappa,r1) 
1252:           if(r2 .gt. DHswitch2)then 
1253:            RD1=r1-DHswitch 
1254:            f_over_r=f_over_r-C2T*(2*COEF2*RD1+3*COEF3*RD1**2) 
1255:            energy=energy-COEF2*RD1**2+COEF3*RD1**3 
1256:           elseif(r2 .lt. DHswitch2)then 
1257:          ! normal DH term 
1258:           else 
1259:          ! things should have fallen in one of the previous two... 
1260:            write(*,*) 'something went wrong with DH switching function' 
1261:           endif 
1262:           f_over_r=f_over_r*sqrt(rm2) 
1263:           grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx 
1264:           grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy 
1265:           grad(C1*3)   = grad(C1*3)   + f_over_r * dz 
1266:           grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx 
1267:           grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy 
1268:           grad(C2*3)   = grad(C2*3)   - f_over_r * dz 
1269:          ENDIF 
1270:        ENDDO 
1271: !$OMP END DO 
1272: !$OMP END PARALLEL1223: !$OMP END PARALLEL
1273:  
1274:  
1275:  
1276:       END1224:       END
1277: 1225: 
1278: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^1226: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^
1279: 1227: 
1280: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1228: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1281: !* SBMPR computes the forces due to position restraints               *1229: !* SBMPR computes the forces due to position restraints               *
1282: !**********************************************************************1230: !**********************************************************************
1283: 1231: 
1284:       subroutine SBMPR(x,y,z,grad, energy, natoms, 1232:       subroutine SBMPR(x,y,z,grad, energy, natoms, 
1285:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)1233:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0