hdiff output

r30685/SBM.f 2016-07-06 15:35:55.279491277 +0100 r30684/SBM.f 2016-07-06 15:35:55.619495874 +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=30000) 11:         parameter(NSBMMAX=25000)
 12:       INTEGER MAXSEP,MAXSEPSYS 12:       INTEGER MAXSEP,MAXSEPSYS
 13:       PARAMETER (MAXSEP=20) 13:       PARAMETER (MAXSEP=20)
 14:       LOGICAL NNCINC 14:       LOGICAL NNCINC
 15:       DIMENSION NNCINC(NSBMMAX*MAXSEP) 15:       DIMENSION NNCINC(NSBMMAX*MAXSEP)
 16:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 16:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX),
 17:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF(NSBMMAX*5*6), 17:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF(NSBMMAX*5*6),
 18:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 18:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
 19:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2), 19:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2),
 20:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3), 20:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3),
 21:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5), 21:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5),
 26:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS), 26:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS),
 27:      Q NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC) 27:      Q NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC)
 28:  28: 
 29:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3) 29:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3)
 30:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX) 30:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX)
 31:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 31:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 32:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 32:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
 33:      Q SBMPRK,SBMPRX,CONCOEF 33:      Q SBMPRK,SBMPRX,CONCOEF
 34:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC, 34:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC,
 35:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1, 35:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,
 36:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,MAXSEPSYS 36:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC
 37:         common /logical/ NNCINC 37:         common /logical/ NNCINC
 38:  38: 
 39:         if(NATOMS.gt. NSBMMAX)then 39:         if(NATOMS.gt. NSBMMAX)then
 40:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 40:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 41:         STOP 41:         STOP
 42:         endif 42:         endif
 43:  43: 
 44:        if(.NOT.CALLED)then 44:        if(.NOT.CALLED)then
 45:         NEXCLUDE=MAXEXCLUSIONS 45:         NEXCLUDE=MAXEXCLUSIONS
 46:         NEXCLUDEELEC=MAXEXCLUSIONS 46:         NEXCLUDEELEC=MAXEXCLUSIONS
 49:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 49:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 50:         write(*,*)  '   For more information on SMOG models, see:' 50:         write(*,*)  '   For more information on SMOG models, see:'
 51:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 51:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 52:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 52:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 53:         write(*,*) 53:         write(*,*)
 54:          54:         
 55:  55: 
 56:        call SBMinit(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP, 56:        call SBMinit(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,
 57:      Q PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE, 57:      Q PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE,
 58:      Q NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC, 58:      Q NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,
 59:      Q SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,CONTACTTYPE, 59:      Q SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,NSBMMAX,CONTACTTYPE,
 60:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 60:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 61:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or. 61:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or.
 62:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then 62:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then
 63:         write(*,*) 'ERROR: arrays not large large enough for the system' 63:         write(*,*) 'increase array size'
 64:         write(*,*) '   entries=',NBA, NTA,NPA,NC  64:         write(*,*) NBA, NTA,NPA,NC 
 65:         stop 65:         stop
 66:         endif  66:         endif 
 67:         CALLED=.TRUE. 67:         CALLED=.TRUE.
 68:         endIF 68:         endIF
 69:  69: 
 70:  70: 
 71: ! call the energy routine 71: ! call the energy routine
 72:       call calc_energy_SBM(qo,natoms,GRAD,energy,Ib1,Ib2,Rb,bK,IT,JT,KT, 72:       call calc_energy_SBM(qo,natoms,NSBMMAX,GRAD,energy,Ib1,Ib2,Rb,bK,IT,JT,KT,
 73:      Q ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF,  73:      Q ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF, 
 74:      Q NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC, 74:      Q NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,
 75:      Q MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,STT,NCcut, 75:      Q MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,STT,NCcut,
 76:      Q CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK, 76:      Q CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,
 77:      Q SBMPRX) 77:      Q SBMPRX)
 78:       IF (STEST) THEN 78:       IF (STEST) THEN
 79:          PRINT '(A)','ERROR - second derivatives not available' 79:          PRINT '(A)','ERROR - second derivatives not available'
 80:          STOP 80:          STOP
 81:       ENDIF 81:       ENDIF
 82:       return 82:       return
 86: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 86: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 87: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 87: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 88: !* startt then the velocities are assigned, otherwise, they are read   * 88: !* startt then the velocities are assigned, otherwise, they are read   *
 89: !* by selecting 2, or generated by selecting 3                         * 89: !* by selecting 2, or generated by selecting 3                         *
 90: !*********************************************************************** 90: !***********************************************************************
 91:  91: 
 92:       subroutine SBMINIT(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk, 92:       subroutine SBMINIT(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,
 93:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF,  93:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF, 
 94:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNEXL1ELEC,NNEXL2ELEC,MAXEXCLUSIONSELEC, 94:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNEXL1ELEC,NNEXL2ELEC,MAXEXCLUSIONSELEC,
 95:      Q NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON, 95:      Q NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,
 96:      Q NCswitch,NCcut,STT,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch, 96:      Q NCswitch,NCcut,STT,NSBMMAX,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,
 97:      Q DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 97:      Q DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 98:       USE KEY 98:       USE KEY
 99:       USE COMMONS, only: ATMASS 99:       USE COMMONS, only: ATMASS
100:       implicit NONE100:       implicit NONE
101: 101: 
102:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,102:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,
103:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 103:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 
104:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,104:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,NSBMMAX,
105:      Q MAXEXCLUSIONS,MAXEXCLUSIONSNEW,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,MAXSEPSYS105:      Q MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,MAXSEPSYS
106: 106: 
107:       DOUBLE PRECISION :: Rb(2*NATOMS), bK(2*NATOMS), ANTC(2*NATOMS), 107:       DOUBLE PRECISION :: Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 
108:      Q Tk(2*NATOMS), PK(3*NATOMS), PHISBM(3*NATOMS),108:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX),
109:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS)109:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS)
110: !      DOUBLE PRECISION CONCOEF110: !      DOUBLE PRECISION CONCOEF
111:       DOUBLE PRECISION CONCOEF (NC*6)111:       DOUBLE PRECISION CONCOEF (NC*6)
112:       INTEGER  Ib1(2*NATOMS), Ib2(2*NATOMS), IT(2*NATOMS), JT(2*NATOMS), 112:       INTEGER  Ib1(2*NSBMMAX), Ib2(2*NSBMMAX), IT(2*NSBMMAX), JT(2*NSBMMAX), 
113:      Q KT(2*NATOMS),IP(3*NATOMS),JP(3*NATOMS),KP(3*NATOMS),113:      Q KT(2*NSBMMAX),IP(3*NSBMMAX),JP(3*NSBMMAX),KP(3*NSBMMAX),
114:      Q LP(3*NATOMS),IC(NATOMS*5),JC(NATOMS*5), 114:      Q LP(3*NSBMMAX),IC(NSBMMAX*5),JC(NSBMMAX*5), 
115:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE115:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE
116:       INTEGER PHITYPE(3*NATOMS),M116:       INTEGER PHITYPE(3*NSBMMAX),M
117:       integer AA,BB,ANTEMP117:       integer AA,BB,ANTEMP
118:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)118:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)
119:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)119:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)
120:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)120:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)
121:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)121:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)
122:       DOUBLE PRECISION dx,dy,dz122:       DOUBLE PRECISION dx,dy,dz
123:       double precision PI,ARG1,ARG2,ARG3,ARG4123:       double precision PI
124:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut124:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
125:       integer CONTACTTYPE(NATOMS*5),TMPARG125:       integer CONTACTTYPE(NSBMMAX*5),TMPARG
126:   
127:       character TMPCHAR126:       character TMPCHAR
128:       integer TMPINT,NUMCHARGES,nexc,I1,I2127:       integer TMPINT,NUMCHARGES,nexc,I1,I2
129:       double precision TMPREAL,concentration128:       double precision TMPREAL,concentration
130:       LOGICAL NNCINC(NATOMS*MAXSEP),TARR(NATOMS)129:       LOGICAL NNCINC(NATOMS*MAXSEP)
131: 130: 
132:       INTEGER SBMPRN, SBMPRI(NATOMS)131:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)
133:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)132:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)
134:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS133:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS
135:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)134:       DIMENSION EXCLUSIONS(NATOMS,MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)
136: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 135: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
137: !$OMP PARALLEL136: !$OMP PARALLEL
138: !$    NTHREADS = OMP_GET_NUM_THREADS()137: !$    NTHREADS = OMP_GET_NUM_THREADS()
139: !$      TID = OMP_GET_THREAD_NUM()138: !$      TID = OMP_GET_THREAD_NUM()
140: !$    if(TID .eq. 0)then 139: !$    if(TID .eq. 0)then 
141: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS140: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS
142: !$    endif141: !$    endif
143: !$OMP END PARALLEL142: !$OMP END PARALLEL
144: 143: 
145:       DO I=1,NATOMS144:       DO I=1,NATOMS
146:         DO J=1,MAXSEP145:         DO J=1,MAXSEP
147:           NNCINC((I-1)*MAXSEP+J)=.TRUE.146:           NNCINC(I*MAXSEP-J+1)=.TRUE.
148:         enddo147:         enddo
149:       enddo148:       enddo
150: 149: 
151:       pi = 3.14159265358979323846264338327950288419716939937510150:       pi = 3.14159265358979323846264338327950288419716939937510
152:       MaxCon=NATOMS*5151:       MaxCon=NATOMS*5
153:       do i=1,NATOMS152:       do i=1,NATOMS
154:         NUMOFEXCLUSIONS(I)=0153:         NUMOFEXCLUSIONS(I)=0
155:       enddo154:       enddo
156: 155: 
157: ! These lines read in the parameters.156: ! These lines read in the parameters.
186:         185:         
187:         read(30,*) 186:         read(30,*) 
188:         read(30,*) NC187:         read(30,*) NC
189:         write(*,*) NC, 'contacts'        188:         write(*,*) NC, 'contacts'        
190:           if(NC .gt. MaxCon)then189:           if(NC .gt. MaxCon)then
191:              write(*,*) 'too many contacts'190:              write(*,*) 'too many contacts'
192:              STOP191:              STOP
193:           endif192:           endif
194: 193: 
195:         do i=1, NC194:         do i=1, NC
196:           read(30, *) IC(i), JC(i), CONTACTTYPE(I),ARG1, ARG2195:           read(30, *) IC(i), JC(i), CONTACTTYPE(I),Sigma, EpsC
197:           !! Add more contact types196:           !! Add more contact types
198:           TMPARG=0197:           TMPARG=0
199:           IF(CONTACTTYPE(I) .NE. 5)THEN198:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,TMPARG)
200:             CALL INCLUDEEXCLUSIONS(NATOMS,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,TMPARG)199:           if(CONTACTTYPE(I) .eq. 1)then
201:           ENDIF200:             CONCOEF(I*6-5)=EPSC*SIGMA**12
202:           if(CONTACTTYPE(I) .eq. 1)then ! 6-12 interaction201:             CONCOEF(I*6-4)=2*EPSC*SIGMA**6
203:             CONCOEF(I*6-5)=ARG2*ARG1**12202:           elseif(CONTACTTYPE(I) .eq. 2)then
204:             CONCOEF(I*6-4)=2*ARG2*ARG1**6203:             CONCOEF(I*6-5)=5*EPSC*SIGMA**12
205:           elseif(CONTACTTYPE(I) .eq. 2)then ! 10-12204:             CONCOEF(I*6-4)=6*EPSC*SIGMA**10
206:             CONCOEF(I*6-5)=5*ARG2*ARG1**12 
207:             CONCOEF(I*6-4)=6*ARG2*ARG1**10 
208:           elseif(CONTACTTYPE(I) .eq. 5)then ! bare Gaussian 
209:             CONCOEF(I*6-5)=ARG1 
210:             CONCOEF(I*6-4)=ARG2 
211:             read(30, *) ARG1 ! need one additional parameter 
212:             CONCOEF(I*6-3)=1/2*ARG1**2 
213:           elseif(CONTACTTYPE(I) .eq. 6)then ! Gaussian with wall 
214:             CONCOEF(I*6-5)=ARG1 
215:             CONCOEF(I*6-4)=ARG2 
216:             read(30, *) ARG1,ARG2 ! need two additional parameters 
217:             CONCOEF(I*6-3)=1/2*ARG1**2 
218:             CONCOEF(I*6-2)=ARG2/CONCOEF(I*6-5) 
219:           elseif(CONTACTTYPE(I) .eq. 7)then ! Dual gaussian 
220:             CONCOEF(I*6-5)=ARG1 
221:             CONCOEF(I*6-4)=ARG2 
222:             read(30, *) ARG1,ARG2,ARG3,ARG4 ! need four additional parameters 
223:             CONCOEF(I*6-3)=1/2*ARG1**2 
224:             CONCOEF(I*6-2)=ARG2 
225:             CONCOEF(I*6-1)=2*ARG3**2 
226:             CONCOEF(I*6)=ARG4/CONCOEF(I*6-5) 
227:           else205:           else
228:             write(*,*) 'ERROR: Unrecognized contacttype:',CONTACTTYPE(I)206:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.'
229:             STOP 
230:           endif207:           endif
231:         end do208:         end do
232: ! make a function to check and add209: ! make a function to check and add
233: 210: 
234:           read(30,*)211:           read(30,*)
235:           read(30,*) nBA212:           read(30,*) nBA
236:           write(*,*) nBA, 'bonds'213:           write(*,*) nBA, 'bonds'
237:         do i=1, nBA214:         do i=1, nBA
238:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)215:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)
239:           !!! ADD BONDS OF TYPE 6216:           !!! ADD BONDS OF TYPE 6
240:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP .AND. bondtype .EQ. 1)THEN217:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP)THEN
241:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS GREATER THAN MAXSEP'218:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS LESS THAN MAXSEP'
242:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'219:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
243:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'220:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
244:           ENDIF221:           ENDIF
245:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN222:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN
246:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          223:             CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
247:           ENDIF224:           ENDIF
248:         end do225:         end do
249: 226: 
250:           read(30,*)227:           read(30,*)
251:           read(30,*) nTA228:           read(30,*) nTA
252:           write(*,*) nTA, 'angles'229:           write(*,*) nTA, 'angles'
253:         do i=1, nTA230:         do i=1, nTA
254:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)231:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)
255:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          232:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
256:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          233:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
257:           CALL INCLUDEEXCLUSIONS(NATOMS,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          234:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
258:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.235:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.
259:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.236:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.
260:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN237:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN
261:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS LESS THAN MAXSEP'238:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS LESS THAN MAXSEP'
262:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'239:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
263:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'240:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
264:           ENDIF241:           ENDIF
265:         ENDDO242:         ENDDO
266: 243: 
267:           read(30,*) 244:           read(30,*) 
268:           read(30,*) nPA245:           read(30,*) nPA
269:           write(*,*) nPA, 'dihedrals'246:           write(*,*) nPA, 'dihedrals'
270: 247: 
271: ! this reads in the dihedral angles and calculates the cosines and sines248: ! this reads in the dihedral angles and calculates the cosines and sines
272: ! in order to make the force and energy calculations easier, later.249: ! in order to make the force and energy calculations easier, later.
273:         do i=1, npA250:         do i=1, npA
274:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)251:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)
275:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),JP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          252:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IP(I),JP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
276:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          253:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
277:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          254:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
278:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          255:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
279:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          256:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
280:           CALL INCLUDEEXCLUSIONS(NATOMS,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          257:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
281: 258: 
282:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.259:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.
283:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.260:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.
284:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.261:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.
285:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.262:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.
286:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.263:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.
287:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN264:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN
288:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS LESS THAN MAXSEP'265:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS LESS THAN MAXSEP'
289:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'266:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
290:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'267:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
291:           ENDIF268:           ENDIF
292:         ENDDO 269:         ENDDO 
293: 270: 
294:         read(30,*)271:         read(30,*)
295:         read(30,*) nexc272:         read(30,*) nexc
296:         write(*,*) nexc, 'exculusions'273:         write(*,*) nexc, 'exculusions'
297: 274: 
298:       do i=1, nexc275:       do i=1, nexc
299:         read(30,*) I1, I2276:         read(30,*) I1, I2
300:         CALL INCLUDEEXCLUSIONS(NATOMS,I1,I2,EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,TMPARG)          277:         CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,I1,I2,EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
301:       enddo278:       enddo
302: 279: 
303: ! now that we know what to exclude, let's make a static list of pairs to280: ! now that we know what to exclude, let's make a static list of pairs to
304: ! exclude281: ! exclude
 282: ! At this point, we will change the meaning of maxexclusions, which is
 283: ! actually NEXCLUDE
305:         NT=MAXEXCLUSIONS*NATOMS      284:         NT=MAXEXCLUSIONS*NATOMS      
306:         NTE=MAXEXCLUSIONSELEC*NATOMS      285:         NTE=MAXEXCLUSIONSELEC*NATOMS      
307:         MAXEXCLUSIONSNEW=0  286:         MAXEXCLUSIONS=0  
308:         MAXEXCLUSIONSELEC=0  287:         MAXEXCLUSIONSELEC=0  
309:         DO I=1,NATOMS288:         DO I=1,NATOMS
310:           DO J=1,NUMOFEXCLUSIONS(I)289:           DO J=1,NUMOFEXCLUSIONS(I)
311:             MAXEXCLUSIONSNEW=MAXEXCLUSIONSNEW+1290:             MAXEXCLUSIONS=MAXEXCLUSIONS+1
312:             IF(MAXEXCLUSIONSNEW .GT. NT)THEN291:             IF(MAXEXCLUSIONS .GT. NT)THEN
313:               write(*,*) 'ERROR: TOO MANY EXCLUSIONS, IN TOTAL.'292:               write(*,*) 'ERROR: TOO MANY EXCLUSIONS, IN TOTAL.'
314:             ENDIF293:             ENDIF
315:             M=EXCLUSIONS((I-1)*MAXEXCLUSIONS+J)294:             M=EXCLUSIONS(I,J)
316:             NNEXL1(MAXEXCLUSIONSNEW)=I295:             NNEXL1(MAXEXCLUSIONS)=I
317:             NNEXL2(MAXEXCLUSIONSNEW)=M296:             NNEXL2(MAXEXCLUSIONS)=M
 297:           
318:             IF(SBMCHARGE(I) .NE. 0 .AND. SBMCHARGE(M) .NE. 0)THEN298:             IF(SBMCHARGE(I) .NE. 0 .AND. SBMCHARGE(M) .NE. 0)THEN
319:               MAXEXCLUSIONSELEC=MAXEXCLUSIONSELEC+1299:               MAXEXCLUSIONSELEC=MAXEXCLUSIONSELEC+1
320:               IF(MAXEXCLUSIONSELEC .GT. NTE)THEN300:               IF(MAXEXCLUSIONSELEC .GT. NTE)THEN
321:                 write(*,*) 'ERROR: TOO MANY EXCLUSIONS (ELEC), IN TOTAL.'301:                 write(*,*) 'ERROR: TOO MANY EXCLUSIONS (ELEC), IN TOTAL.'
322:               ENDIF302:               ENDIF
323:               NNEXL1ELEC(MAXEXCLUSIONSELEC)=I303:               NNEXL1ELEC(MAXEXCLUSIONSELEC)=I
324:               NNEXL2ELEC(MAXEXCLUSIONSELEC)=M304:               NNEXL2ELEC(MAXEXCLUSIONSELEC)=M
325:             ENDIF305:             ENDIF
326:           ENDDO306:           ENDDO
327:         ENDDO307:         ENDDO
328:         STT=reps*rsig**12308:         STT=reps*rsig**12
329: ! At this point, we will change the meaning of maxexclusions, which is309: 
330: ! actually NEXCLUDE 
331:         MAXEXCLUSIONS=MAXEXCLUSIONSNEW 
332:         do I=1,NATOMS310:         do I=1,NATOMS
333:           TARR(I)=.FALSE.311:           TARR(I)=0
334:         enddo312:         enddo
335: 313: 
336: ! read in position restraints314: ! read in position restraints
337:        read(30,*) 315:        read(30,*) 
338:        read(30,*) SBMPRN316:        read(30,*) SBMPRN
339:        write(*,*) SBMPRN, 'position restraints'317:        write(*,*) SBMPRN, 'position restraints'
340:        do I=1,SBMPRN318:        do I=1,SBMPRN
341:             J=(I-1)*6319:             J=(I-1)*6
342:             K=J/2320:             K=J/2
343:             read(30,*) SBMPRI(I),SBMPRK(J+1),SBMPRK(J+2),SBMPRK(J+3),321:             read(30,*) SBMPRI(I),SBMPRK(J+1),SBMPRK(J+2),SBMPRK(J+3),
344:      Q         SBMPRK(J+4),SBMPRK(J+5),SBMPRK(J+6),322:      Q         SBMPRK(J+4),SBMPRK(J+5),SBMPRK(J+6),
345:      Q         SBMPRX(K+1),SBMPRX(K+2),SBMPRX(K+3)323:      Q         SBMPRX(K+1),SBMPRX(K+2),SBMPRX(K+3)
346:             if(TARR(SBMPRI(I)))then324:             if(TARR(SBMPRI(I)) .ne. 0)then
347:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)325:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
348:                 STOP326:                 STOP
349:             endif327:             endif
350:             TARR(SBMPRI(I))=.TRUE.328:             TARR(SBMPRI(I))=1
351:         if(SBMPRK(J+4) .ne. 0 .or. SBMPRK(J+5) .ne. 0 .or. SBMPRK(J+6) .ne.0)then  329:         if(SBMPRK(J+4) .ne. 0 .or. SBMPRK(J+5) .ne. 0 .or. SBMPRK(J+6) .ne.0)then  
352:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'330:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'
353:         endif331:         endif
354:        enddo 332:        enddo 
355:        close(30)333:        close(30)
356:       END 334:       END 
357: 335: 
358: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^336: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
359: 337: 
360: 338: 
361: !**************************************************************************339: !**************************************************************************
362: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in340: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in
363: ! the model341: ! the model
364: !**************************************************************************342: !**************************************************************************
365: 343: 
366:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,344:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,
367:      Q NNCINC,MAXSEP,MAXSEPSYS)345:      Q NNCINC,MAXSEP,MAXSEPSYS)
368:       IMPLICIT NONE346:       IMPLICIT NONE
369:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,347:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,
370:      Q MAXSEPSYS348:      Q MAXSEPSYS,NSBMMAX
371:       DIMENSION EXCLIST(NATOMS*MAXEXCLUDE),NEXCLUDE(NATOMS)349:       DIMENSION EXCLIST(NATOMS,MAXEXCLUDE),NEXCLUDE(NATOMS)
372:       LOGICAL  NNCINC350:       LOGICAL  NNCINC
373:       DIMENSION NNCINC(NATOMS*MAXSEP)351:       DIMENSION NNCINC(NATOMS*MAXSEP)
374:       INTEGER I,N,M,AA,BB,AB352:       INTEGER I,N,M,AA,BB,AB
375:       LOGICAL INCLUDED353:       LOGICAL INCLUDED
376: 354: 
377: ! If the atoms are within MAXSEP of each other, then keep track via355: ! If the atoms are within MAXSEP of each other, then keep track via
378: ! NNCINC356: ! NNCINC
379:       AA=min(ATOM1,ATOM2)357:       AA=min(ATOM1,ATOM2)
380:       BB=max(ATOM1,ATOM2)358:       BB=max(ATOM1,ATOM2)
381:       AB=BB-AA359:       AB=BB-AA
382:       IF(AB .le. MAXSEP)THEN360:       IF(AB .le. MAXSEP)THEN
383:         NNCINC((AA-1)*MAXSEP+AB)=.FALSE.361:         NNCINC((AA-1)*MAXSEP+AB)=.FALSE.
384:         IF(AB .gt. MAXSEPSYS)THEN362:         IF(AB .gt. MAXSEPSYS)THEN
385:           MAXSEPSYS=AB363:           MAXSEPSYS=AB
386:         ENDIF364:         ENDIF
387:       ELSE365:       ELSE
388:         INCLUDED = .FALSE.366:         INCLUDED = .FALSE.
389:         N=MAX(ATOM1,ATOM2)367:         N=MAX(ATOM1,ATOM2)
390:         M=MIN(ATOM1,ATOM2)368:         M=MIN(ATOM1,ATOM2)
391:         DO I =1,NEXCLUDE(M)369:         DO I =1,NEXCLUDE(M)
392:           IF(EXCLIST((M-1)*MAXEXCLUDE+I) .eq. N)THEN370:           IF(EXCLIST(M,I) .eq. N)THEN
393:             INCLUDED = .TRUE.371:             INCLUDED = .TRUE.
394:           ENDIF372:           ENDIF
395:         ENDDO 373:         ENDDO 
396:   374:   
397:         if(.NOT. INCLUDED)THEN375:         if(.NOT. INCLUDED)THEN
398:           if(NEXCLUDE(M) .EQ. MAXEXCLUDE)THEN376:           if(NEXCLUDE(M) .EQ. MAXEXCLUDE)THEN
399:             write(*,*) 'ERROR: TOO MANY EXCLUSIONS WITH ATOM',M377:             write(*,*) 'ERROR: TOO MANY EXCLUSIONS WITH ATOM',M
400:           ENDIF378:           ENDIF
401:           NEXCLUDE(M)=NEXCLUDE(M)+1379:           NEXCLUDE(M)=NEXCLUDE(M)+1
402:           EXCLIST((M-1)*MAXEXCLUDE+NEXCLUDE(M))=N380:           EXCLIST(M,NEXCLUDE(M))=N
403:         ENDIF381:         ENDIF
404:       ENDIF382:       ENDIF
405:       END383:       END
406: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^384: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
407: 385: 
408: 386: 
409: C387: C
410: C Calculate the Forces and energies388: C Calculate the Forces and energies
411: C389: C
412:       subroutine CALC_ENERGY_SBM(qo,natoms,GRAD,energy,Ib1,Ib2,390:       subroutine CALC_ENERGY_SBM(qo,natoms,NSBMMAX,GRAD,energy,Ib1,Ib2,
413:      Q Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,391:      Q Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,
414:      Q CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,NNEXL2ELEC,392:      Q CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,NNEXL2ELEC,
415:      Q NEXCLUSIONSELEC,NCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,393:      Q NEXCLUSIONSELEC,NCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,
416:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,394:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
417:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)395:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
418: 396: 
419:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,397:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,
420:      Q NEXCLUSIONSELEC398:      Q NSBMMAX,NEXCLUSIONSELEC
421: 399: 
422:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY400:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
423:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)401:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
424:       LOGICAL NNCINC402:       LOGICAL NNCINC
425:       DIMENSION NNCINC(NATOMS*MAXSEP)403:       DIMENSION NNCINC(NATOMS*MAXSEP)
426: 404: 
427:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 405:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
428:      Q PHISBM(NPA),NCswitch,NCcut,STT, 406:      Q PHISBM(NPA),NCswitch,NCcut,STT, 
429:      Q SBMCHARGE(NATOMS)407:      Q SBMCHARGE(NATOMS)
430:         DOUBLE PRECISION  CONCOEF(NC*6)408:         DOUBLE PRECISION  CONCOEF(NC*6)
446:         grad(j+1) = 0.0424:         grad(j+1) = 0.0
447:         grad(j+2) = 0.0425:         grad(j+2) = 0.0
448:         grad(j+3) = 0.0426:         grad(j+3) = 0.0
449:       enddo427:       enddo
450: 428: 
451:       energy = 0.0429:       energy = 0.0
452:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)430:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
453:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)431:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
454:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,432:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
455:      Q PHITYPE,PHISBM,NPA)433:      Q PHITYPE,PHISBM,NPA)
456:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC,CONCOEF,434:       call SBMContacts(x,y,z,grad, energy,natoms,NSBMMAX,IC,JC,CONCOEF,
457:      Q NC,CONTACTTYPE)435:      Q NC,CONTACTTYPE)
458:       call SBMNonContacts(x,y,z,grad,energy,natoms,NNEXL1,NNEXL2,436:       call SBMNonContacts(x,y,z,grad,energy,natoms,NSBMMAX,NNEXL1,NNEXL2,
459:      Q NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)437:      Q NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)
460:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,438:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,
461:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)439:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
462:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)440:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
463: 441: 
464:       end442:       end
465: 443: 
466: 444: 
467: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<445: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
468: !* SBMBonds  computes the hookean force and energy between chosen atoms *446: !* SBMBonds  computes the hookean force and energy between chosen atoms *
756: 734: 
757: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^735: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^
758: 736: 
759: 737: 
760: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<738: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
761: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *739: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *
762: !* 10-12 or 6-12 potential                                              *740: !* 10-12 or 6-12 potential                                              *
763: !***********************************************************************741: !***********************************************************************
764: 742: 
765:       subroutine SBMcontacts(x,y,z,grad,energy,743:       subroutine SBMcontacts(x,y,z,grad,energy,
766:      Q NATOMS,IC,JC,CONCOEF,NC,CONTACTTYPE)744:      Q NATOMS,NSBMMAX,IC,JC,CONCOEF,NC,CONTACTTYPE)
767:       USE KEY745:       USE KEY
768:       implicit NONE746:       implicit NONE
769:       integer I, NATOMS,NC747:       integer I, NATOMS,NC,NSBMMAX
770: 748: 
771:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 749:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
772:      Q grad(3*NATOMS), energy750:      Q grad(3*NATOMS), energy
773:       DOUBLE PRECISION dx,dy,dz,CO1,CO2,CO3,CO4,CO5,CO6,E1,G1,G2,E1D,G1D,G2D751:       DOUBLE PRECISION dx,dy,dz,CO1,CO2
774: 752: 
775:       integer C1, C2 753:       integer C1, C2 
776:       DOUBLE PRECISION  r2,rm2,rm10,rm12,rm6,f_over_r,r1 754:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 
777: 755: 
778:         DOUBLE PRECISION CONCOEF(NC*6)756:         DOUBLE PRECISION CONCOEF(NC*6)
779:         INTEGER IC(NC), JC(NC),CONTACTTYPE(NC)757:         INTEGER IC(NC), JC(NC),CONTACTTYPE(NC)
780: 758: 
781: ! type 1 is 6-12759: ! type 1 is 6-12
782: ! type 2 is 12-12760: ! type 2 is 12-12
783: ! implement type 3 as gaussian761: ! implement type 3 as gaussian
784: ! implement type 4 as dual-basin gaussian762: ! implement type 4 as dual-basin gaussian
785: ! question: Should type be part of the interaction, or should they be organized763: ! question: Should type be part of the interaction, or should they be organized
786: ! into separate lists?764: ! into separate lists?
787: ! type 1 is 6-12 interaction765: ! type 1 is 6-12 interaction
788: !$OMP PARALLEL PRIVATE(I,CO1,CO2,CO3,CO4,CO5,CO6,C1,C2,DX,DY,DZ,766: !$OMP PARALLEL PRIVATE(I,CO1,CO2,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R) REDUCTION(+:ENERGY,grad)
789: !$OMP& R2,RM2,RM6,RM10,RM12,R1,E1,G1,G2,E1D,G1D,G2D,F_OVER_R) REDUCTION(+:ENERGY,grad) 
790: !$OMP DO767: !$OMP DO
791:           do i=1, NC768:           do i=1, NC
792:             C1 = IC(i)769:             C1 = IC(i)
793:             C2 = JC(i)770:             C2 = JC(i)
794:             dx = X(C1) - X(C2)771:             dx = X(C1) - X(C2)
795:             dy = Y(C1) - Y(C2)772:             dy = Y(C1) - Y(C2)
796:             dz = Z(C1) - Z(C2)773:             dz = Z(C1) - Z(C2)
797:             r2 = dx**2 + dy**2 + dz**2774:             r2 = dx**2 + dy**2 + dz**2
 775:             rm2 = 1.0/r2
 776: 
798:             if(CONTACTTYPE(I) .eq. 1)then777:             if(CONTACTTYPE(I) .eq. 1)then
799:               ! type 1 is 6-12 interaction778:               ! type 1 is 6-12 interaction
800:               rm2 = 1.0/r2 
801:               rm6 = rm2**3779:               rm6 = rm2**3
802:               CO1=CONCOEF(I*6-5)780:               CO1=CONCOEF(I*6-5)
803:               CO2=CONCOEF(I*6-4)781:               CO2=CONCOEF(I*6-4)
804:               ENERGY = ENERGY + rm6*(CO1*rm6-CO2)782:               ENERGY = ENERGY + rm6*(CO1*rm6-CO2)
805:               f_over_r = -12.0*rm6*(CO1*rm6-CO2/2.0)*rm2783:               f_over_r = -12.0*rm6*(CO1*rm6-CO2/2.0)*rm2
806:             elseif(CONTACTTYPE(I) .eq. 2)then784:             elseif(CONTACTTYPE(I) .eq. 2)then
807:               ! type 2 is 10-12 interaction785:               ! type 2 is 10-12 interaction
808:               rm2 = 1.0/r2 
809:               rm10 = rm2**5786:               rm10 = rm2**5
810:               CO1=CONCOEF(I*6-5)787:               CO1=CONCOEF(I*6-5)
811:               CO2=CONCOEF(I*6-4)788:               CO2=CONCOEF(I*6-4)
812:               ENERGY = ENERGY + rm10*(CO1*rm2-CO2)789:               ENERGY = ENERGY + rm10*(CO1*rm2-CO2)
813:               f_over_r = -60.0*rm10*(CO1*rm2/5.0-CO2/6.0)*rm2790:               f_over_r = -60.0*rm10*(CO1*rm2/5.0-CO2/6.0)*rm2
814:             elseif(CONTACTTYPE(I) .eq. 5)then 
815:               ! type 5 is a simple gaussian interaction (no excluded 
816:               ! volume) 
817:               r1 = sqrt(r2) 
818:               CO1=CONCOEF(I*6-5) !A 
819:               CO2=CONCOEF(I*6-4) !mu 
820:               CO3=CONCOEF(I*6-3) !1/(2*sigma**2) 
821:               E1= -CO1*exp(-CO3*(r1-CO2)**2) 
822:               ENERGY = ENERGY + E1 
823:               f_over_r = -E1*2*CO3*(r1-CO2)/r1 
824:             elseif(CONTACTTYPE(I) .eq. 6)then 
825:               ! type 6 is a gaussian interaction with excluded 
826:               ! volume included 
827:               r1 = sqrt(r2) 
828:               rm2 = 1.0/r2 
829:               rm12 = rm2**6 
830:               CO1=CONCOEF(I*6-5) !A 
831:               CO2=CONCOEF(I*6-4) !mu 
832:               CO3=CONCOEF(I*6-3) !1/(2*sigma**2) 
833:               CO4=CONCOEF(I*6-2) !a/A 
834:               E1 =CO1*(1+CO4*rm12) 
835:               G1 =1-exp(-CO3*(r1-CO2)**2) 
836:               E1D=-12*CO1*CO4*rm12*rm2 
837:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2) 
838:               ENERGY = ENERGY + (E1*G1-CO1) 
839:               f_over_r = E1*G1D + G1*E1D  
840:             elseif(CONTACTTYPE(I) .eq. 7)then 
841:               ! type 7 is a double gaussian interaction with excluded 
842:               ! volume included 
843:               rm2 = 1.0/r2 
844:               rm12 = rm2**6 
845:               r1 = sqrt(r2) 
846:               CO1=CONCOEF(I*6-5) !A 
847:               CO2=CONCOEF(I*6-4) !mu1 
848:               CO3=CONCOEF(I*6-3) !1/(2*sigma1**2) 
849:               CO4=CONCOEF(I*6-2) !mu2 
850:               CO5=CONCOEF(I*6-1) !1/(2*sigma2**2) 
851:               CO6=CONCOEF(I*6)   ! a/A (excluded volume amplitude/gaussian amplitude) 
852:               E1=1+CO6*rm12 
853:               G1=1-exp(-CO3*(r1-CO2)**2) 
854:               G2=1-exp(-CO5*(r1-CO4)**2) 
855:               E1D= -12*CO6*rm12*rm2 
856:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2) 
857:               G2D=CO5*2*(1-CO4/r1)*exp(-CO5*(r1-CO4)**2) 
858:               ENERGY = ENERGY + CO1*( E1 * G1 * G2)  
859:               f_over_r = CO1*(E1D*G1*G2  
860:      Q        +  E1 * G1D * G2 
861:      Q        +  E1 * G1 * G2D) 
862:             else791:             else
863:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)792:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)
864:               STOP793:               STOP
865:             endif794:             endif
866: 795: 
867:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx796:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
868:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy797:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
869:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz798:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz
870:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx799:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
871:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy800:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
877: 806: 
878:       end807:       end
879: 808: 
880: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^809: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
881: 810: 
882: 811: 
883: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<812: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
884: !* SBMNonContacts computes the forces due to non native contacts       *813: !* SBMNonContacts computes the forces due to non native contacts       *
885: !**********************************************************************814: !**********************************************************************
886: 815: 
887:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,NNEXL1,NNEXL2,816:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,NSBMMAX,NNEXL1,NNEXL2,
888:      Q  NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)817:      Q  NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)
889:         USE KEY818:         USE KEY
890:         implicit NONE819:         implicit NONE
891:         integer I, N, J, AN, NATOMS820:         integer I, N, J, AN, NATOMS
892: 821: 
893: 822: 
894:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),823:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
895:      Q  grad(3*NATOMS), energy,STT,SA,SB,SC,ET824:      Q  grad(3*NATOMS), energy,STT,SA,SB,SC
896: 825: 
897:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj826:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj
898:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 827:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
899:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,MAXSEPSYS828:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,NSBMMAX,MAXSEPSYS
900:         LOGICAL NNCINC829:         LOGICAL NNCINC
901:         DIMENSION NNCINC(NATOMS*MAXSEP)830:         DIMENSION NNCINC(NATOMS*MAXSEP)
902: !$    INTEGER  OMP_GET_NUM_THREADS,831: !$    INTEGER  OMP_GET_NUM_THREADS,
903: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS832: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
904: 833: 
905:         INTEGER NNEXL1(NEXCLUSIONS)834:         INTEGER NNEXL1(NEXCLUSIONS)
906:         INTEGER NNEXL2(NEXCLUSIONS)835:         INTEGER NNEXL2(NEXCLUSIONS)
907:         DOUBLE PRECISION dx,dy,dz836:         DOUBLE PRECISION dx,dy,dz
908:         integer tempN, alpha837:         integer tempN, alpha
909:         double precision Rdiff,Vfunc,Ffunc838:         double precision Rdiff,Vfunc,Ffunc
926:         GRIDSIZE=NCcut*1.01855:         GRIDSIZE=NCcut*1.01
927:         Rcut2=NCcut**2856:         Rcut2=NCcut**2
928:         Rswitch2=NCswitch**2857:         Rswitch2=NCswitch**2
929: 858: 
930:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 859:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
931:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )860:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
932: 861: 
933:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)862:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)
934: 863: 
935:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)864:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)
 865: 
936: ! grid the system866: ! grid the system
937:         minX=10000000867:         minX=10000000
938:         minY=10000000868:         minY=10000000
939:         minZ=10000000869:         minZ=10000000
940:         870:         
941:         maxX=-10000000871:         maxX=-10000000
942:         maxY=-10000000872:         maxY=-10000000
943:         maxZ=-10000000873:         maxZ=-10000000
944: 874: 
945:         do i=1,NATOMS875:         do i=1,NATOMS
997:                 endif927:                 endif
998:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i928:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
999:         enddo929:         enddo
1000: 930: 
1001: 931: 
1002:            tempN=0932:            tempN=0
1003: 933: 
1004: ! add a second loop that goes over charged atoms934: ! add a second loop that goes over charged atoms
1005: 935: 
1006: !$OMP PARALLEL936: !$OMP PARALLEL
1007: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,937: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID)  
1008: !$OMP& Xgrid,Ygrid,Zgrid,TID,DINDEX)   
1009: !$OMP& REDUCTION(+:ENERGY,grad)938: !$OMP& REDUCTION(+:ENERGY,grad)
1010:         TID=0939:         TID=0
1011:         NTHREADS=1;940:         NTHREADS=1;
1012: !$      TID = OMP_GET_THREAD_NUM()941: !$      TID = OMP_GET_THREAD_NUM()
1013: !$      NTHREADS = OMP_GET_NUM_THREADS()942: !$      NTHREADS = OMP_GET_NUM_THREADS()
1014:         do i=1+TID,NATOMS,NTHREADS943:         do i=1+TID,NATOMS,NTHREADS
1015: 944: 
1016:           Xgrid=int((X(I)-minX)/gridsize)+1945:           Xgrid=int((X(I)-minX)/gridsize)+1
1017:           Ygrid=int((Y(I)-minY)/gridsize)+1946:           Ygrid=int((Y(I)-minY)/gridsize)+1
1018:           Zgrid=int((Z(I)-minZ)/gridsize)+1947:           Zgrid=int((Z(I)-minZ)/gridsize)+1
1019:           do ii=XGRID-1,XGRID+1948:           do ii=XGRID-1,XGRID+1
1020:            do jj=YGRID-1,YGRID+1949:            do jj=YGRID-1,YGRID+1
1021:             do kk=ZGRID-1,ZGRID+1950:             do kk=ZGRID-1,ZGRID+1
1022:              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.
1023:      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
1024:               do jjj=1,Ngrid(ii,jj,kk)953:               do jjj=1,Ngrid(ii,jj,kk)
 954: 
1025:                C2=grid(ii,jj,kk,jjj)955:                C2=grid(ii,jj,kk,jjj)
1026:                DINDEX=C2-I956:                DINDEX=C2-I
1027:                if(DINDEX .GT. 0)THEN957:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNEXL(C1,C12))) )then
1028:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX)  ) )then958:            !if(DINDEX .gt. MAXSEP .or. (DINDEX .gt. 0 .and. NNEXL(I,DINDEX)))then
1029:                 if(I .EQ. 5178 .AND. C2 .EQ. 5212)THEN959:                if(DINDEX .gt. 0 .AND. ( DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX) ) ) )then
1030:                 endif960:                 dx = X(I) - X(C2)
1031:                  dx = X(I) - X(C2)961:                 dy = Y(I) - Y(C2)
1032:                  dy = Y(I) - Y(C2)962:                 dz = Z(I) - Z(C2)
1033:                  dz = Z(I) - Z(C2)963:                 r2 = dx**2 + dy**2 + dz**2
1034:                  r2 = dx**2 + dy**2 + dz**2964:                 if(r2 .le. Rcut2)then
1035:                  if(r2 .le. Rcut2)then965:                  rm2 = 1/r2
1036:                   rm2 = 1/r2966:                  rm14 = rm2**7
1037:                   rm14 = rm2**7967:                  energy=energy+STT*rm2**6+SC
1038:                   energy=energy+STT*rm2**6+SC968:                  f_over_r=-STT*12.0*rm14
1039:                   f_over_r=-STT*12.0*rm14969:                  RD1=sqrt(r2)-NCswitch
1040:                   RD1=sqrt(r2)-NCswitch970:                  if(r2 .gt. Rswitch2)then
1041:                   if(r2 .gt. Rswitch2)then971:                   f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
1042:                    f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)972:                   energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4
1043:                    energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4973:                  elseif(r2 .lt. Rswitch2)then
1044:                   elseif(r2 .lt. Rswitch2)then974:                 ! normal repulsive term
1045:                  ! normal repulsive term975:                  else
1046:                   else976:                 ! things should have fallen in one of the previous two...
1047:                  ! things should have fallen in one of the previous two...977:                   write(*,*) 'something went wrong with switching function'
1048:                    write(*,*) 'something went wrong with switching function'978:                  endif
1049:                   endif979:                  grad(I*3-2) = grad(I*3-2) + f_over_r * dx
1050:                   grad(I*3-2) = grad(I*3-2) + f_over_r * dx980:                  grad(I*3-1) = grad(I*3-1) + f_over_r * dy
1051:                   grad(I*3-1) = grad(I*3-1) + f_over_r * dy981:                  grad(I*3)   = grad(I*3)   + f_over_r * dz
1052:                   grad(I*3)   = grad(I*3)   + f_over_r * dz982:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
1053:                   grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx983:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
1054:                   grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy984:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz
1055:                   grad(C2*3)   = grad(C2*3)   - f_over_r * dz985:                  endif
1056:                   endif986:                endif
1057:                 endif 
1058:                endif 
1059:               enddo987:               enddo
1060:              endif988:              endif
1061:             enddo989:             enddo
1062:            enddo990:            enddo
1063:           enddo991:           enddo
1064:          enddo992:          enddo
1065: 993: 
1066: !$OMP DO994: !$OMP DO
1067:        DO I=1,NEXCLUSIONS995:        DO I=1,NEXCLUSIONS
1068:          C1=NNEXL1(I) 996:          C1=NNEXL1(I) 
1069:          C2=NNEXL2(I)997:          C2=NNEXL2(I)
1070:          dx = X(C1) - X(C2)998:          dx = X(C1) - X(C2)
1071:          dy = Y(C1) - Y(C2)999:          dy = Y(C1) - Y(C2)
1072:          dz = Z(C1) - Z(C2)1000:          dz = Z(C1) - Z(C2)
1073:          r2 = dx**2 + dy**2 + dz**21001:          r2 = dx**2 + dy**2 + dz**2
1074:          if(r2 .le. Rcut2)then1002:          if(r2 .le. Rcut2)then
1075:           rm2 = 1/r21003:           rm2 = 1/r2
1076:           rm14 = rm2**71004:           rm14 = rm2**7
1077:            energy=energy-STT*rm2**6-SC1005:           energy=energy-STT*rm2**6-SC
1078:            ET=-STT*rm2**6-SC1006:           f_over_r=-STT*12.0*rm14
1079:            f_over_r=-STT*12.0*rm141007:           RD1=sqrt(r2)-NCswitch
1080:            RD1=sqrt(r2)-NCswitch 
1081:           if(r2 .gt. Rswitch2)then1008:           if(r2 .gt. Rswitch2)then
1082:            f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)1009:            f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
1083:            energy=energy-SA/3.0*RD1**3-SB/4.0*RD1**41010:            energy=energy-SA/3.0*RD1**3-SB/4.0*RD1**4
1084:            ET=ET-SA/3.0*RD1**3-SB/4.0*RD1**41011:           elseif(r2 .lt. Rswitch2)then
1085: !          elseif(r2 .lt. Rswitch2)then 
1086:          ! normal repulsive term1012:          ! normal repulsive term
1087: !          else1013:           else
1088:          ! things should have fallen in one of the previous two...1014:          ! things should have fallen in one of the previous two...
1089: !           write(*,*) 'something went wrong with switching function'1015:            write(*,*) 'something went wrong with switching function'
1090:           endif1016:           endif
1091:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx1017:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx
1092:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy1018:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy
1093:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz1019:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz
1094:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx1020:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx
1095:           grad(C2*3-1) = grad(C2*3-1) + f_over_r * dy1021:           grad(C2*3-1) = grad(C2*3-1) + f_over_r * dy
1096:           grad(C2*3)   = grad(C2*3)   + f_over_r * dz1022:           grad(C2*3)   = grad(C2*3)   + f_over_r * dz
1097:          endif1023:          endif
1098:       ENDDO1024:       ENDDO
1099: !$OMP END DO1025: !$OMP END DO
1156:         ! dimensions of grid1082:         ! dimensions of grid
1157:         parameter (maxgrid=15)1083:         parameter (maxgrid=15)
1158:         dimension Ngrid(maxgrid,maxgrid,maxgrid),1084:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
1159:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)1085:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
1160:         integer MaxGridX,MaxGridY,MaxGridZ1086:         integer MaxGridX,MaxGridY,MaxGridZ
1161:         double precision gridsize,RD11087:         double precision gridsize,RD1
1162:         double precision minX,minY,minZ,maxX,maxY,maxZ1088:         double precision minX,minY,minZ,maxX,maxY,maxZ
1163:         integer Xgrid,Ygrid,Zgrid,TID1089:         integer Xgrid,Ygrid,Zgrid,TID
1164:         double precision C2T1090:         double precision C2T
1165: 1091: 
 1092: 
1166:       if(SBMCHARGEON(1) .eq. 1) RETURN1093:       if(SBMCHARGEON(1) .eq. 1) RETURN
1167: 1094: 
1168:         diff=DHcut-DHswitch1095:         diff=DHcut-DHswitch
1169:         alpha=121096:         alpha=12
1170: 1097: 
1171: 1098: 
1172:         GRIDSIZE=DHCUT*1.011099:         GRIDSIZE=DHCUT*1.01
1173:         DHcut2=DHcut**21100:         DHcut2=DHcut**2
1174:         DHswitch2=DHswitch**21101:         DHswitch2=DHswitch**2
1175: 1102: 
1386:            DX=X(J)-SBMPRX(I2+1)1313:            DX=X(J)-SBMPRX(I2+1)
1387:            DY=Y(J)-SBMPRX(I2+2)1314:            DY=Y(J)-SBMPRX(I2+2)
1388:            DZ=Z(J)-SBMPRX(I2+3)1315:            DZ=Z(J)-SBMPRX(I2+3)
1389:            K1=SBMPRK(I1+1)1316:            K1=SBMPRK(I1+1)
1390:            K2=SBMPRK(I1+2)1317:            K2=SBMPRK(I1+2)
1391:            K3=SBMPRK(I1+3)1318:            K3=SBMPRK(I1+3)
1392:            K4=SBMPRK(I1+4)1319:            K4=SBMPRK(I1+4)
1393:            K5=SBMPRK(I1+5)1320:            K5=SBMPRK(I1+5)
1394:            K6=SBMPRK(I1+6)1321:            K6=SBMPRK(I1+6)
1395:            1322:            
 1323: 
1396:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +1324:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +
1397:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ1325:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ
1398: 1326: 
1399:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ 1327:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ 
1400:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ1328:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ
1401:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY1329:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY
1402:         enddo1330:         enddo
1403: !$OMP END DO1331: !$OMP END DO
1404: !$OMP END PARALLEL1332: !$OMP END PARALLEL
1405: 1333: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0