hdiff output

r30696/SBM.f 2016-07-06 15:36:00.107556559 +0100 r30695/SBM.f 2016-07-06 15:36:00.447561158 +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:       USE SBMDATA  
  4:       implicit NONE  3:       implicit NONE
  5:       INTEGER NATOMS,i,j  4:       INTEGER NATOMS,i,j
  6:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)
  7:       DOUBLE PRECISION ENERGY  6:       DOUBLE PRECISION ENERGY
  8:   7: 
  9:       LOGICAL :: CALLED=.FALSE.  8:       LOGICAL :: CALLED=.FALSE.
 10:       LOGICAL GTEST, STEST  9:       LOGICAL GTEST, STEST
  10:         integer NSBMMAX
  11:         parameter(NSBMMAX=12000)
  12:       INTEGER MAXSEP,MAXSEPSYS
  13:       PARAMETER (MAXSEP=40)
  14:       INTEGER MAXATOMTYPES
  15:       INTEGER ATOMTYPES
  16:       PARAMETER (MAXATOMTYPES=10)
  17:       DOUBLE PRECISION STT,SA, SB, SC
  18:       DIMENSION STT(MAXATOMTYPES*MAXATOMTYPES),ATOMTYPES(NSBMMAX)
  19:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
  20:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
  21:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
  22:       INTEGER NNCINC
  23:       DIMENSION NNCINC(NSBMMAX*MAXSEP)
  24:       DOUBLE PRECISION  Rb(5*NSBMMAX),bK(5*NSBMMAX),ANTC(2*NSBMMAX),
  25:      Q Tk(2*NSBMMAX),PK(3*NSBMMAX),PHISBM(3*NSBMMAX),CONCOEF(NSBMMAX*4*6),
  26:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
  27:       INTEGER  Ib1(NSBMMAX*5), Ib2(NSBMMAX*5), IT(NSBMMAX*2), JT(NSBMMAX*2),
  28:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3),
  29:      Q LP(NSBMMAX*3), IC(NSBMMAX*4), JC(NSBMMAX*4),
  30:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX)
  31:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NEXCLUDE,NEXCLUDEELEC
  32:       PARAMETER (MAXEXCLUSIONS=20)
  33:       PARAMETER (MAXEXCLUSIONSELEC=5)
  34:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),
  35:      Q        NNEXL2(NSBMMAX*MAXEXCLUSIONS),
  36:      Q        NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),
  37:      Q        NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC)
  38: 
  39:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3)
  40:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX)
  41:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
  42:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
  43:      Q SBMPRK,SBMPRX,CONCOEF,SA,SB,SC
  44:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC,
  45:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,
  46:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,MAXSEPSYS,
  47:      Q ATOMTYPES
 11:  48: 
 12:         if(.NOT.CALLED)then 49:         if(NATOMS.gt. NSBMMAX)then
  50:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
  51:         STOP
  52:         endif
  53: 
  54:        if(.NOT.CALLED)then
  55:         NEXCLUDE=MAXEXCLUSIONS
  56:         NEXCLUDEELEC=MAXEXCLUSIONS
  57:         MAXSEPSYS=0
  58:       DO I=1,NATOMS
  59:         DO J=1,MAXSEP
  60:           NNCINC((I-1)*MAXSEP+J)=0
  61:         enddo
  62:       enddo
 13:  63: 
 14:         write(*,*) 64:         write(*,*)
 15:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 65:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 16:         write(*,*)  '   For more information on SMOG models, see:' 66:         write(*,*)  '   For more information on SMOG models, see:'
 17:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 67:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 18:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 68:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 19:         write(*,*) 69:         write(*,*)
  70:         
 20:  71: 
 21:         call SBMinit(NATOMS) 72:        call SBMinit(NATOMS,ATOMTYPES,MAXATOMTYPES,Ib1,Ib2,Rb,bK,IT,JT,
  73:      Q KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1,
  74:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,
  75:      Q MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,
  76:      Q STT,SA,SB,SC,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,
  77:      Q SBMPRI,SBMPRK,SBMPRX)
  78:         if(NBA .gt. NSBMMAX*5 .or. NTA .gt. NSBMMAX*2 .or.
  79:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 4*NSBMMAX)then
  80:         write(*,*) 'ERROR: arrays not large large enough for the system'
  81:         write(*,*) '   entries=',NBA, NTA,NPA,NC 
  82:         stop
  83:         endif 
 22:         CALLED=.TRUE. 84:         CALLED=.TRUE.
 23:         endIF 85:         endIF
 24:  86: 
 25:  87: 
 26: ! call the energy routine 88: ! call the energy routine
 27:       call calc_energy_SBM(qo,natoms,GRAD,energy) 89:       call calc_energy_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,energy,
  90:      Q Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,
  91:      Q JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,
  92:      Q NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,
  93:      Q SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,
  94:      Q DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
  95: 
 28:  96: 
 29:       IF (STEST) THEN 97:       IF (STEST) THEN
 30:          PRINT '(A)','ERROR - second derivatives not available' 98:          PRINT '(A)','ERROR - second derivatives not available'
 31:          STOP 99:          STOP
 32:       ENDIF100:       ENDIF
 33:       return101:       return
 34:       end102:       end
 35: 103: 
 104: 
 36: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<105: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 37: !* SBMinit() reads the atom positions from file.  If 1 is selected for *106: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 38: !* startt then the velocities are assigned, otherwise, they are read   *107: !* startt then the velocities are assigned, otherwise, they are read   *
 39: !* by selecting 2, or generated by selecting 3                         *108: !* by selecting 2, or generated by selecting 3                         *
 40: !***********************************************************************109: !***********************************************************************
 41: 110: 
 42:       subroutine SBMINIT(NATOMS)111:       subroutine SBMINIT(NATOMS,ATOMTYPES,MAXATOMTYPES,Ib1,Ib2,Rb,bK,IT,
 112:      Q JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF, 
 113:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNEXL1ELEC,NNEXL2ELEC,MAXEXCLUSIONSELEC,
 114:      Q NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,
 115:      Q NCswitch,NCcut,STT,SA,SB,SC,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,
 116:      Q DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 43:       USE KEY117:       USE KEY
 44:       USE SBMDATA 
 45:       USE COMMONS, only: ATMASS118:       USE COMMONS, only: ATMASS
 46:       USE GENRIGID, only: RIGIDINIT 
 47:       implicit NONE119:       implicit NONE
 48: 120: 
 49:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,121:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,
 50:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 122:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 
 51:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,123:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,
 52:      Q ATYPE,POS124:      Q MAXEXCLUSIONS,MAXEXCLUSIONSNEW,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,
 53:       DOUBLE PRECISION :: Sigma, EpsC125:      Q MAXSEPSYS,ATOMTYPES(NATOMS),ATYPE,POS,MAXATOMTYPES,NATOMTYPES
 54:       DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:) :: TMPSTT126:       DOUBLE PRECISION :: Rb(2*NATOMS), bK(2*NATOMS), ANTC(2*NATOMS), 
 127:      Q Tk(2*NATOMS), PK(3*NATOMS), PHISBM(3*NATOMS),
 128:      Q Sigma, EpsC, NCswitch,NCcut,TMPSTT(MAXATOMTYPES),
 129:      Q STT(MAXATOMTYPES**2),SBMCHARGE(NATOMS)
 55:       DOUBLE PRECISION RDIFF,ALPHA130:       DOUBLE PRECISION RDIFF,ALPHA
 56:       INTEGER BONDTYPE131:       DOUBLE PRECISION CONCOEF(NC*6)
 57:       INTEGER M132:       INTEGER  Ib1(5*NATOMS), Ib2(5*NATOMS), IT(2*NATOMS), JT(2*NATOMS), 
 133:      Q KT(2*NATOMS),IP(3*NATOMS),JP(3*NATOMS),KP(3*NATOMS),
 134:      Q LP(3*NATOMS),IC(NATOMS*4),JC(NATOMS*4), 
 135:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE
 136:       INTEGER PHITYPE(3*NATOMS),M
 58:       integer AA,BB,ANTEMP137:       integer AA,BB,ANTEMP
 138:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)
 139:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)
 140:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)
 141:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)
 59:       DOUBLE PRECISION dx,dy,dz142:       DOUBLE PRECISION dx,dy,dz
 60:       double precision PI,ARG1,ARG2,ARG3,ARG4143:       double precision PI,ARG1,ARG2,ARG3,ARG4
 61:       DOUBLE PRECISION RSig, Reps,DC144:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
 62:       integer TMPARG145:       integer CONTACTTYPE(NATOMS*5),TMPARG
 146:       DOUBLE PRECISION SA, SB, SC
 147:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
 148:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
 149:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
 150: 
 63:       character TMPCHAR151:       character TMPCHAR
 64:       integer TMPINT,nexc,I1,I2152:       integer TMPINT,NUMCHARGES,nexc,I1,I2
 65:       double precision TMPREAL,concentration153:       double precision TMPREAL,concentration
 66:       LOGICAL TARR(NATOMS)154:       LOGICAL TARR(NATOMS)
 155:       INTEGER, DIMENSION(NATOMS*MAXSEP) ::  NNCINC
 156:       INTEGER SBMPRN, SBMPRI(NATOMS)
 157:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)
 67:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS158:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS
 68:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC 
 69:       PARAMETER(MAXEXCLUSIONS=20) 
 70:       PARAMETER(MAXEXCLUSIONSELEC=5) 
 71:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)159:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)
 72:       INTEGER IOSTATUS,J1,J2,NRIGIDBODY,ATOMID 
 73:       CHARACTER CHECK1 
 74: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 160: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
 75: !$OMP PARALLEL161: !$OMP PARALLEL
 76: !$    NTHREADS = OMP_GET_NUM_THREADS()162: !$    NTHREADS = OMP_GET_NUM_THREADS()
 77: !$      TID = OMP_GET_THREAD_NUM()163: !$      TID = OMP_GET_THREAD_NUM()
 78: !$    if(TID .eq. 0)then 164: !$    if(TID .eq. 0)then 
 79: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS165: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS
 80: !$    endif166: !$    endif
 81: !$OMP END PARALLEL167: !$OMP END PARALLEL
 82: 168: 
 83:       ALLOCATE(XSBM(NATOMS)) 
 84:       ALLOCATE(YSBM(NATOMS)) 
 85:       ALLOCATE(ZSBM(NATOMS)) 
 86:       ALLOCATE(ATOMTYPES(NATOMS)) 
 87: 169: 
 88:       pi = 3.14159265358979323846264338327950288419716939937510170:       pi = 3.14159265358979323846264338327950288419716939937510
 89:       MaxCon=NATOMS*5171:       MaxCon=NATOMS*5
 90:       do i=1,NATOMS172:       do i=1,NATOMS
 91:         NUMOFEXCLUSIONS(I)=0173:         NUMOFEXCLUSIONS(I)=0
 92:       enddo174:       enddo
 93: 175: 
 94:  
 95:         MAXSEP=40 
 96:         MAXSEPSYS=0 
 97:       ALLOCATE(NNCINC(NATOMS*MAXSEP)) 
 98:       DO I=1,NATOMS 
 99:         DO J=1,MAXSEP 
100:           NNCINC((I-1)*MAXSEP+J)=0 
101:         enddo 
102:       enddo 
103:  
104: ! These lines read in the parameters.176: ! These lines read in the parameters.
105:         open(30, file='SBM.INP', status='old', access='sequential')177:         open(30, file='SBM.INP', status='old', access='sequential')
106:         write(*,*) 178:         write(*,*) 
107:         write(*,*) 'Reading the forcefield from SBM.INP'179:         write(*,*) 'Reading the forcefield from SBM.INP'
108:         write(*,*) 180:         write(*,*) 
109:         read(30,*)181:         read(30,*)
110:         read(30,*)182:         read(30,*)
111:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut183:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut
112:         read(30,*)184:         read(30,*)
113:         read(30,*) NCswitch, NCcut185:         read(30,*) NCswitch, NCcut
114:         write(*,*) 'NCswitch,NCcut'186:         write(*,*) 'NCswitch,NCcut'
115:         write(*,'(2F10.5)') NCswitch, NCcut187:         write(*,'(2F10.5)') NCswitch, NCcut
116:         write(*,*) 'Reading electrostatic parameters'188:         write(*,*) 'Reading electrostatic parameters'
117:         write(*,'(5F10.5)') PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut189:         write(*,'(5F10.5)') PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut
118:         read(30,*) NATOMTYPES190:         read(30,*) NATOMTYPES
119:         ALLOCATE(TMPSTT(NATOMTYPES)) 
120:         do I=1,NATOMTYPES191:         do I=1,NATOMTYPES
121:           read(30,*) TMPARG,rsig,reps192:           read(30,*) TMPARG,rsig,reps
122:           if(TMPARG .ne. I)then193:           if(TMPARG .ne. I)then
123:             write(*,*) 'ERROR: atomtypes must be sequential integers'194:             write(*,*) 'ERROR: atomtypes must be sequential integers'
124:             STOP  195:             STOP  
125:           endif196:           endif
126:           TMPSTT(TMPARG)=reps*rsig**12197:           TMPSTT(TMPARG)=reps*rsig**12
127:         enddo198:         enddo
128:         ! using combination rule 1199:         ! using combination rule 1
129: 200: 
130:         Rdiff=NCcut-NCswitch201:         Rdiff=NCcut-NCswitch
131:         alpha=12202:         alpha=12
132:         ALLOCATE(STT(NATOMTYPES**2)) 203: 
133:         ALLOCATE(SA(NATOMTYPES**2))  
134:         ALLOCATE(SB(NATOMTYPES**2))  
135:         ALLOCATE(SC(NATOMTYPES**2))  
136:         do I=1,NATOMTYPES204:         do I=1,NATOMTYPES
137:           do J=1,NATOMTYPES205:           do J=1,NATOMTYPES
138:             POS=(I-1)*NATOMTYPES+J206:             POS=(I-1)*MAXATOMTYPES+J
139:             STT(POS)=sqrt(TMPSTT(I)*TMPSTT(J))207:             STT(POS)=sqrt(TMPSTT(I)*TMPSTT(J))
140:             SB(POS)=-1.0/Rdiff**3*( 2*alpha*STT(POS)/NCcut**(alpha+1)  + 208:             SB(POS)=-1.0/Rdiff**3*( 2*alpha*STT(POS)/NCcut**(alpha+1)  + 
141:      Q         (alpha)*(alpha+1)*STT(POS)*Rdiff/NCcut**(alpha+2))209:      Q         (alpha)*(alpha+1)*STT(POS)*Rdiff/NCcut**(alpha+2))
142:             SA(POS)=-(alpha*(alpha+1)*STT(POS)/NCcut**(alpha+2)+210:             SA(POS)=-(alpha*(alpha+1)*STT(POS)/NCcut**(alpha+2)+
143:      Q         3*SB(POS)*Rdiff**2)/(2*Rdiff)211:      Q         3*SB(POS)*Rdiff**2)/(2*Rdiff)
144:             SC(POS)=-(STT(POS)/NCcut**alpha +212:             SC(POS)=-(STT(POS)/NCcut**alpha +
145:      Q         SA(POS)/3.0*Rdiff**3+SB(POS)/4.0*Rdiff**4)213:      Q         SA(POS)/3.0*Rdiff**3+SB(POS)/4.0*Rdiff**4)
146:           enddo 214:           enddo 
147:         enddo215:         enddo
148: 216: 
149:         DEALLOCATE(TMPSTT) 
150: 217: 
151: 218: 
152:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is in units219:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is in units
153:         ! A^-1220:         ! A^-1
154:         KAPPA=0.32*sqrt(CONCENTRATION)221:         KAPPA=0.32*sqrt(CONCENTRATION)
155:         PREFACTOR=PREFACTOR/DC222:         PREFACTOR=PREFACTOR/DC
156:         read(30,*) ANtemp223:         read(30,*) ANtemp
157:         write(*,*) ANtemp, 'atoms'224:         write(*,*) ANtemp, 'atoms'
158:         if(NATOMS .EQ. ANTEMP)THEN225:         NUMCHARGES=1
159:           write(*,*) 'Number of atoms in SBM.INP and odata are not consistent' 
160:         ENDIF 
161:         NUMOFSBMCHARGES=0 
162:         ALLOCATE(SBMCHARGE(ANTEMP)) 
163:         ALLOCATE(SBMCHARGEON(ANTEMP)) 
164:         do i=1, ANtemp226:         do i=1, ANtemp
165:           read(30,*) TMPINT,ATYPE,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)227:           read(30,*) TMPINT,ATYPE,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)
166:             if(ATYPE .gt. NATOMTYPES)THEN228:             if(ATYPE .gt. NATOMTYPES)THEN
167:               write(*,*) 'ERROR: Unknown atomtype',ATYPE229:               write(*,*) 'ERROR: Unknown atomtype',ATYPE
168:             endif230:             endif
169:             ATOMTYPES(I)=ATYPE231:             ATOMTYPES(I)=ATYPE
170:           if(SBMCHARGE(i) .ne. 0)then232:           if(SBMCHARGE(i) .ne. 0)then
171:              NUMOFSBMCHARGES=NUMOFSBMCHARGES+1233:              NUMCHARGES=NUMCHARGES+1
172:              SBMCHARGEON(NUMOFSBMCHARGES)=i234:              SBMCHARGEON(NUMCHARGES)=i
173:           endif 235:           endif 
174:         end do236:         end do
 237:         SBMCHARGEON(1)=NUMCHARGES
175:         238:         
176:         read(30,*) 239:         read(30,*) 
177:         read(30,*) NC240:         read(30,*) NC
178:         write(*,*) NC, 'contacts'        241:         write(*,*) NC, 'contacts'        
179:           if(NC .gt. MaxCon)then242:           if(NC .gt. MaxCon)then
180:              write(*,*) 'too many contacts'243:              write(*,*) 'too many contacts'
181:              STOP244:              STOP
182:           endif245:           endif
183: 246: 
184:         ALLOCATE(IC(NC)) 
185:         ALLOCATE(JC(NC)) 
186:         ALLOCATE(CONTACTTYPE(NC)) 
187:         ALLOCATE(CONCOEF(NC*6)) 
188:  
189:         do i=1, NC247:         do i=1, NC
190:           read(30, *) IC(i), JC(i), CONTACTTYPE(I),ARG1, ARG2248:           read(30, *) IC(i), JC(i), CONTACTTYPE(I),ARG1, ARG2
191:           !! Add more contact types249:           !! Add more contact types
192:           TMPARG=0250:           TMPARG=0
193:           IF(CONTACTTYPE(I) .NE. 5)THEN251:           IF(CONTACTTYPE(I) .NE. 5)THEN
194:             CALL INCLUDEEXCLUSIONS(NATOMS,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.FALSE.)252:             CALL INCLUDEEXCLUSIONS(NATOMS,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,TMPARG)
195:           ENDIF253:           ENDIF
196:           if(CONTACTTYPE(I) .eq. 1)then ! 6-12 interaction254:           if(CONTACTTYPE(I) .eq. 1)then ! 6-12 interaction
197:             CONCOEF(I*6-5)=ARG2*ARG1**12255:             CONCOEF(I*6-5)=ARG2*ARG1**12
198:             CONCOEF(I*6-4)=2*ARG2*ARG1**6256:             CONCOEF(I*6-4)=2*ARG2*ARG1**6
199:           elseif(CONTACTTYPE(I) .eq. 2)then ! 10-12257:           elseif(CONTACTTYPE(I) .eq. 2)then ! 10-12
200:             CONCOEF(I*6-5)=5*ARG2*ARG1**12258:             CONCOEF(I*6-5)=5*ARG2*ARG1**12
201:             CONCOEF(I*6-4)=6*ARG2*ARG1**10259:             CONCOEF(I*6-4)=6*ARG2*ARG1**10
202:           elseif(CONTACTTYPE(I) .eq. 5)then ! bare Gaussian260:           elseif(CONTACTTYPE(I) .eq. 5)then ! bare Gaussian
203:             CONCOEF(I*6-5)=ARG1261:             CONCOEF(I*6-5)=ARG1
204:             CONCOEF(I*6-4)=ARG2262:             CONCOEF(I*6-4)=ARG2
219:             CONCOEF(I*6-1)=1/(2*ARG3**2)277:             CONCOEF(I*6-1)=1/(2*ARG3**2)
220:             CONCOEF(I*6)=ARG4/CONCOEF(I*6-5)278:             CONCOEF(I*6)=ARG4/CONCOEF(I*6-5)
221:           else279:           else
222:             write(*,*) 'ERROR: Unrecognized contacttype:',CONTACTTYPE(I)280:             write(*,*) 'ERROR: Unrecognized contacttype:',CONTACTTYPE(I)
223:             STOP281:             STOP
224:           endif282:           endif
225:         end do283:         end do
226: ! make a function to check and add284: ! make a function to check and add
227: 285: 
228:           read(30,*)286:           read(30,*)
229:           read(30,*) NBA287:           read(30,*) nBA
230:           write(*,*) NBA, 'bonds'288:           write(*,*) nBA, 'bonds'
231:           ALLOCATE(IB1(NBA)) 
232:           ALLOCATE(IB2(NBA)) 
233:           ALLOCATE(RB(NBA)) 
234:           ALLOCATE(BK(NBA)) 
235:         do i=1, nBA289:         do i=1, nBA
236:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)290:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)
237:           !!! ADD BONDS OF TYPE 6291:           !!! ADD BONDS OF TYPE 6
238:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP .AND. bondtype .EQ. 1)THEN292:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP .AND. bondtype .EQ. 1)THEN
239:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS GREATER THAN MAXSEP'293:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS GREATER THAN MAXSEP'
240:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'294:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
241:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'295:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
242:           ENDIF296:           ENDIF
243:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN297:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN
244:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          298:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
245:           ELSEIF(bondtype .eq. 6)THEN 
246:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.FALSE.)           
247:           ENDIF299:           ENDIF
248:         end do300:         end do
249: 301: 
250:           read(30,*)302:           read(30,*)
251:           read(30,*) NTA303:           read(30,*) nTA
252:           write(*,*) NTA, 'angles'304:           write(*,*) nTA, 'angles'
253:         ALLOCATE(IT(NTA)) 
254:         ALLOCATE(JT(NTA)) 
255:         ALLOCATE(KT(NTA)) 
256:         ALLOCATE(ANTC(NTA)) 
257:         ALLOCATE(TK(NTA)) 
258:         do i=1, nTA305:         do i=1, nTA
259:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), TK(i)306:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)
260:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          307:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
261:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          308:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
262:           CALL INCLUDEEXCLUSIONS(NATOMS,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          309:           CALL INCLUDEEXCLUSIONS(NATOMS,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
263:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.310:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.
264:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.311:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.
265:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN312:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN
266:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS GREATER THAN MAXSEP'313:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS GREATER THAN MAXSEP'
267:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'314:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
268:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'315:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
269:           ENDIF316:           ENDIF
270:         ENDDO317:         ENDDO
271: 318: 
272:           read(30,*) 319:           read(30,*) 
273:           read(30,*) NPA320:           read(30,*) nPA
274:           write(*,*) NPA, 'dihedrals'321:           write(*,*) nPA, 'dihedrals'
275:           ALLOCATE(IP(NPA))322: 
276:           ALLOCATE(JP(NPA)) 
277:           ALLOCATE(KP(NPA)) 
278:           ALLOCATE(LP(NPA)) 
279:           ALLOCATE(PK(NPA)) 
280:           ALLOCATE(PHITYPE(NPA)) 
281:           ALLOCATE(PHISBM(NPA)) 
282: ! this reads in the dihedral angles and calculates the cosines and sines323: ! this reads in the dihedral angles and calculates the cosines and sines
283: ! in order to make the force and energy calculations easier, later.324: ! in order to make the force and energy calculations easier, later.
284:         do i=1, npA325:         do i=1, npA
285:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)326:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)
286:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),JP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          327:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),JP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
287:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          328:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
288:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          329:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
289:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          330:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
290:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          331:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
291:           CALL INCLUDEEXCLUSIONS(NATOMS,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.TRUE.)          332:           CALL INCLUDEEXCLUSIONS(NATOMS,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
292: 333: 
293:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.334:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.
294:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.335:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.
295:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.336:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.
296:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.337:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.
297:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.338:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.
298:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN339:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN
299:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS GREATER THAN MAXSEP'340:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS GREATER THAN MAXSEP'
300:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'341:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
301:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'342:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
302:           ENDIF343:           ENDIF
303:         ENDDO 344:         ENDDO 
304: 345: 
305:         read(30,*)346:         read(30,*)
306:         read(30,*) nexc347:         read(30,*) nexc
307:         write(*,*) nexc, 'exculusions'348:         write(*,*) nexc, 'exculusions'
308: 349: 
309:       do i=1, nexc350:       do i=1, nexc
310:         read(30,*) I1, I2351:         read(30,*) I1, I2
311:         CALL INCLUDEEXCLUSIONS(NATOMS,I1,I2,EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,.FALSE.)          352:         CALL INCLUDEEXCLUSIONS(NATOMS,I1,I2,EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,TMPARG)          
312:       enddo353:       enddo
313: 354: 
314: ! now that we know what to exclude, let's make a static list of pairs to355: ! now that we know what to exclude, let's make a static list of pairs to
315: ! exclude356: ! exclude
316:         NEXCLUSIONS=0357:         NT=MAXEXCLUSIONS*NATOMS      
317:         DO I=1,NATOMS358:         NTE=MAXEXCLUSIONSELEC*NATOMS      
318:            NEXCLUSIONS=NEXCLUSIONS+NUMOFEXCLUSIONS(I)359:         MAXEXCLUSIONSNEW=0  
319:         ENDDO360:         MAXEXCLUSIONSELEC=0  
320:         ALLOCATE(NNEXL1(NEXCLUSIONS)) 
321:         ALLOCATE(NNEXL2(NEXCLUSIONS)) 
322:         ALLOCATE(NNEXL1ELEC(NEXCLUSIONS)) 
323:         ALLOCATE(NNEXL2ELEC(NEXCLUSIONS)) 
324:         ! reset counters 
325:         NEXCLUSIONS=0 
326:         NEXCLUSIONSELEC=0 
327:         DO I=1,NATOMS361:         DO I=1,NATOMS
328:           DO J=1,NUMOFEXCLUSIONS(I)362:           DO J=1,NUMOFEXCLUSIONS(I)
329:             NEXCLUSIONS=NEXCLUSIONS+1363:             MAXEXCLUSIONSNEW=MAXEXCLUSIONSNEW+1
 364:             IF(MAXEXCLUSIONSNEW .GT. NT)THEN
 365:               write(*,*) 'ERROR: TOO MANY EXCLUSIONS, IN TOTAL.'
 366:             ENDIF
330:             M=EXCLUSIONS((I-1)*MAXEXCLUSIONS+J)367:             M=EXCLUSIONS((I-1)*MAXEXCLUSIONS+J)
331:             NNEXL1(NEXCLUSIONS)=I368:             NNEXL1(MAXEXCLUSIONSNEW)=I
332:             NNEXL2(NEXCLUSIONS)=M369:             NNEXL2(MAXEXCLUSIONSNEW)=M
333:             IF(SBMCHARGE(I) .NE. 0 .AND. SBMCHARGE(M) .NE. 0)THEN370:             IF(SBMCHARGE(I) .NE. 0 .AND. SBMCHARGE(M) .NE. 0)THEN
334:               NEXCLUSIONSELEC=NEXCLUSIONSELEC+1371:               MAXEXCLUSIONSELEC=MAXEXCLUSIONSELEC+1
335:               NNEXL1ELEC(NEXCLUSIONSELEC)=I372:               IF(MAXEXCLUSIONSELEC .GT. NTE)THEN
336:               NNEXL2ELEC(NEXCLUSIONSELEC)=M373:                 write(*,*) 'ERROR: TOO MANY EXCLUSIONS (ELEC), IN TOTAL.'
 374:               ENDIF
 375:               NNEXL1ELEC(MAXEXCLUSIONSELEC)=I
 376:               NNEXL2ELEC(MAXEXCLUSIONSELEC)=M
337:             ENDIF377:             ENDIF
338:           ENDDO378:           ENDDO
339:         ENDDO379:         ENDDO
340: ! At this point, we will change the meaning of maxexclusions, which is380: ! At this point, we will change the meaning of maxexclusions, which is
341: ! actually NEXCLUDE381: ! actually NEXCLUDE
 382:         MAXEXCLUSIONS=MAXEXCLUSIONSNEW
342:         do I=1,NATOMS383:         do I=1,NATOMS
343:           TARR(I)=.FALSE.384:           TARR(I)=.FALSE.
344:         enddo385:         enddo
345: 386: 
346: ! read in position restraints387: ! read in position restraints
347:        read(30,*) 388:        read(30,*) 
348:        read(30,*) SBMPRN389:        read(30,*) SBMPRN
349:        write(*,*) SBMPRN, 'position restraints'390:        write(*,*) SBMPRN, 'position restraints'
350:  
351:       ALLOCATE(SBMPRI(SBMPRN)) 
352:       ALLOCATE(SBMPRK(SBMPRN*6)) 
353:       ALLOCATE(SBMPRX(SBMPRN*3)) 
354:  
355:        do I=1,SBMPRN391:        do I=1,SBMPRN
356:             J=(I-1)*6392:             J=(I-1)*6
357:             K=J/2393:             K=J/2
358:             read(30,*) SBMPRI(I),SBMPRK(J+1),SBMPRK(J+2),SBMPRK(J+3),394:             read(30,*) SBMPRI(I),SBMPRK(J+1),SBMPRK(J+2),SBMPRK(J+3),
359:      Q         SBMPRK(J+4),SBMPRK(J+5),SBMPRK(J+6),395:      Q         SBMPRK(J+4),SBMPRK(J+5),SBMPRK(J+6),
360:      Q         SBMPRX(K+1),SBMPRX(K+2),SBMPRX(K+3)396:      Q         SBMPRX(K+1),SBMPRX(K+2),SBMPRX(K+3)
361:             if(TARR(SBMPRI(I)))then397:             if(TARR(SBMPRI(I)))then
362:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)398:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
363:                 STOP399:                 STOP
364:             endif400:             endif
365:             TARR(SBMPRI(I))=.TRUE.401:             TARR(SBMPRI(I))=.TRUE.
366:         if(SBMPRK(J+4) .ne. 0 .or. SBMPRK(J+5) .ne. 0 .or. SBMPRK(J+6) .ne.0)then  402:         if(SBMPRK(J+4) .ne. 0 .or. SBMPRK(J+5) .ne. 0 .or. SBMPRK(J+6) .ne.0)then  
367:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'403:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'
368:         endif404:         endif
369:        enddo 405:        enddo 
370:        close(30)406:        close(30)
371:  
372:  
373: ! If RIGIDINIT is on, then we will want to keep information about rigid groups 
374: ! so that we may ignore their interactions later. 
375:  
376:       ALLOCATE (SBMRBG(NATOMS)) 
377:       DO I=1,NATOMS 
378:         SBMRBG(I)=-I 
379:       ENDDO 
380:  
381: ! approach to I/O for rbodyconfig taken from genrigid.f90 
382:       if(RIGIDINIT)THEN 
383:         OPEN(UNIT=222,FILE='rbodyconfig',status='old') 
384:         DO 
385:           READ(222,*,IOSTAT=iostatus) CHECK1 
386:           IF (iostatus<0) THEN 
387:             CLOSE(222) 
388:             EXIT 
389:           ELSE IF (TRIM(ADJUSTL(CHECK1)).EQ.'GROUP') then 
390:             NRIGIDBODY=NRIGIDBODY+1 
391:           ENDIF 
392:         END DO 
393:         CLOSE(222) 
394:     
395:         OPEN(UNIT=222,FILE='rbodyconfig',status='old') 
396:         DO J1 = 1, NRIGIDBODY 
397:           READ(222,*) CHECK1,DUMMY 
398:           DO J2 = 1, DUMMY 
399:             READ(222,*) ATOMID 
400:             SBMRBG(ATOMID)=J1 
401:           ENDDO 
402:         ENDDO 
403:         CLOSE(222) 
404:       ENDIF 
405:  
406:       END 407:       END 
407: 408: 
408: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^409: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
409: 410: 
410: 411: 
411: !**************************************************************************412: !**************************************************************************
412: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in413: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in
413: ! the model414: ! the model
414: !**************************************************************************415: !**************************************************************************
415: 416: 
416:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,417:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,
417:      Q MAXEXCLUDE,CHECKSEP)418:      Q NNCINC,MAXSEP,MAXSEPSYS)
418:       USE SBMDATA 
419:       IMPLICIT NONE419:       IMPLICIT NONE
420:       LOGICAL CHECKSEP420:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,
421:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE421:      Q MAXSEPSYS
422:       DIMENSION EXCLIST(NATOMS*MAXEXCLUDE),NEXCLUDE(NATOMS)422:       DIMENSION EXCLIST(NATOMS*MAXEXCLUDE),NEXCLUDE(NATOMS)
 423:       INTEGER, DIMENSION(NATOMS*MAXSEP) ::  NNCINC
423:       INTEGER I,N,M,AA,BB,AB,POS424:       INTEGER I,N,M,AA,BB,AB,POS
424:       LOGICAL INCLUDED425:       LOGICAL INCLUDED
425: ! If the atoms are within MAXSEP of each other, then keep track via426: ! If the atoms are within MAXSEP of each other, then keep track via
426: ! NNCINC427: ! NNCINC
427:       AA=min(ATOM1,ATOM2)428:       AA=min(ATOM1,ATOM2)
428:       BB=max(ATOM1,ATOM2)429:       BB=max(ATOM1,ATOM2)
429:       AB=BB-AA430:       AB=BB-AA
430:       POS=(AA-1)*MAXSEP+AB431:       POS=(AA-1)*MAXSEP+AB
431:       IF(AB .le. MAXSEP)THEN432:       IF(AB .le. MAXSEP)THEN
432:         NNCINC(POS)=1433:         NNCINC(POS)=1
433:         IF(AB .GT. MAXSEPSYS .AND. CHECKSEP)THEN434:         IF(AB .gt. MAXSEPSYS)THEN
434:           MAXSEPSYS=AB435:           MAXSEPSYS=AB
435:         ENDIF436:         ENDIF
436:       ELSE437:       ELSE
437:         INCLUDED = .FALSE.438:         INCLUDED = .FALSE.
438:         N=MAX(ATOM1,ATOM2)439:         N=MAX(ATOM1,ATOM2)
439:         M=MIN(ATOM1,ATOM2)440:         M=MIN(ATOM1,ATOM2)
440:         DO I =1,NEXCLUDE(M)441:         DO I =1,NEXCLUDE(M)
441:           IF(EXCLIST((M-1)*MAXEXCLUDE+I) .eq. N)THEN442:           IF(EXCLIST((M-1)*MAXEXCLUDE+I) .eq. N)THEN
442:             INCLUDED = .TRUE.443:             INCLUDED = .TRUE.
443:           ENDIF444:           ENDIF
452:         ENDIF453:         ENDIF
453:       ENDIF454:       ENDIF
454: 455: 
455:       END456:       END
456: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^457: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
457: 458: 
458: 459: 
459: C460: C
460: C Calculate the Forces and energies461: C Calculate the Forces and energies
461: C462: C
462:       subroutine CALC_ENERGY_SBM(qo,natoms,GRAD,energy)463:       subroutine CALC_ENERGY_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,
463:       USE SBMDATA464:      Q energy,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,
 465:      Q PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,
 466:      Q NNEXL2ELEC,NEXCLUSIONSELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,
 467:      Q SBMCHARGE,SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,
 468:      Q KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
464: 469: 
465:       INTEGER I,J,NATOMS470:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,
 471:      Q NEXCLUSIONSELEC,ATOMTYPES(NATOMS),MAXATOMTYPES
466: 472: 
467:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY473:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
468: 474:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
 475:       INTEGER NNCINC
 476:       DIMENSION NNCINC(NATOMS*MAXSEP)
 477:       DOUBLE PRECISION SA, SB, SC
 478:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
 479:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
 480:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
 481: 
 482:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
 483:      Q PHISBM(NPA),NCswitch,NCcut,STT(MAXATOMTYPES**2), SBMCHARGE(NATOMS)
 484:         DOUBLE PRECISION  CONCOEF(NC*6)
 485:         INTEGER Ib1(NBA),Ib2(NBA),IT(NTA),JT(NTA),KT(NTA),IP(NPA), 
 486:      Q JP(NPA),KP(NPA),LP(NPA),IC(NC),JC(NC),SBMCHARGEON(NATOMS),
 487:      Q PHITYPE(NPA),CONTACTTYPE(NC)
 488:         INTEGER NNEXL1(NEXCLUSIONS)
 489:         INTEGER NNEXL2(NEXCLUSIONS)
 490:         INTEGER NNEXL1ELEC(NEXCLUSIONSELEC)
 491:         INTEGER NNEXL2ELEC(NEXCLUSIONSELEC)
469:       DOUBLE PRECISION dx,dy,dz492:       DOUBLE PRECISION dx,dy,dz
 493:       INTEGER SBMPRN, SBMPRI(NATOMS)
 494:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)
470:       do i = 1, natoms495:       do i = 1, natoms
471:         j = (i-1)*3496:         j = (i-1)*3
472:         XSBM(i) = qo(j+1)497:         x(i) = qo(j+1)
473:         YSBM(i) = qo(j+2)498:         y(i) = qo(j+2)
474:         ZSBM(i) = qo(j+3)499:         z(i) = qo(j+3)
475:         grad(j+1) = 0.0500:         grad(j+1) = 0.0
476:         grad(j+2) = 0.0501:         grad(j+2) = 0.0
477:         grad(j+3) = 0.0502:         grad(j+3) = 0.0
478:       enddo503:       enddo
479:       energy = 0.0504:       energy = 0.0
480:       call SBMbonds(grad, energy, natoms)505:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
481:       call SBMangl(grad, energy, natoms)506:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
482:       call SBMDihedral(grad, energy, natoms)507:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
483:       call SBMContacts(grad, energy,natoms)508:      Q PHITYPE,PHISBM,NPA)
484:       call SBMNonContacts(grad,energy,natoms)509:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC,CONCOEF,
485:       call SBMDHELEC(grad,energy,natoms)510:      Q NC,CONTACTTYPE)
486:       call SBMPR(grad, energy, natoms)511:       call SBMNonContacts(x,y,z,grad,energy,natoms,ATOMTYPES,
 512:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,
 513:      Q NCswitch,NCcut,STT,SA,SB,SC)
 514:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,
 515:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
 516:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
487:       end517:       end
488: 518: 
489: 519: 
490: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<520: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
491: !* SBMBonds  computes the hookean force and energy between chosen atoms *521: !* SBMBonds  computes the hookean force and energy between chosen atoms *
492: !***********************************************************************522: !***********************************************************************
493: 523: 
494:       subroutine SBMBonds(grad,energy, natoms)524:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)
495:       USE KEY525:       USE KEY
496:       USE SBMDATA 
497:       implicit NONE526:       implicit NONE
498:       integer I2,J2,I23,J23,I,N,J,NATOMS527:       integer I2,J2,outE,I,N,J,NATOMS,NBA,C1
499:       DOUBLE PRECISION grad(3*NATOMS),energy528:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),
 529:      Q energy
500:       DOUBLE PRECISION r2, f, r1530:       DOUBLE PRECISION r2, f, r1
501:       DOUBLE PRECISION dx,dy,dz531:       DOUBLE PRECISION dx,dy,dz
 532:       DOUBLE PRECISION Rb(NBA), bK(NBA)
 533:       INTEGER Ib1(NBA), Ib2(NBA)
502: 534: 
503: !$OMP PARALLEL PRIVATE(I,I2,J2,I23,J23,DX,DY,DZ,R1,R2,F)REDUCTION(+:ENERGY,grad)535: !$OMP PARALLEL PRIVATE(I,I2,J2,DX,DY,DZ,R1,R2,F)REDUCTION(+:ENERGY,grad)
504: !$OMP DO536: !$OMP DO
505:         DO I=1, nBA537:         DO I=1, nBA
506:            I2 = Ib1(I)538:            I2 = Ib1(I)
507:            J2 = Ib2(I)539:            J2 = Ib2(I)
508:            I23=I2*3540:            dx = X(I2) - X(J2)
509:            J23=J2*3541:            dy = Y(I2) - Y(J2)
510:            dx = XSBM(I2) - XSBM(J2)542:            dz = Z(I2) - Z(J2)
511:            dy = YSBM(I2) - YSBM(J2) 
512:            dz = ZSBM(I2) - ZSBM(J2) 
513:            r2 = dx**2 + dy**2 + dz**2543:            r2 = dx**2 + dy**2 + dz**2
514:            r1 = sqrt(r2)544:            r1 = sqrt(r2)
515:            ENERGY = ENERGY + bk(I)*(r1-Rb(I))**2/2.0545:            ENERGY = ENERGY + bk(I)*(r1-Rb(I))**2/2.0
516:            f = -bk(I)*(r1-Rb(I))/r1546:            f = -bk(I)*(r1-Rb(I))/r1
517:            grad(I23-2) = grad(I23-2) - f * dx547:            grad(I2*3-2) = grad(I2*3-2) - f * dx
518:            grad(I23-1) = grad(I23-1) - f * dy548:            grad(I2*3-1) = grad(I2*3-1) - f * dy
519:            grad(I23)   = grad(I23)   - f * dz549:            grad(I2*3)   = grad(I2*3)   - f * dz
520:            grad(J23-2) = grad(J23-2) + f * dx550:            grad(J2*3-2) = grad(J2*3-2) + f * dx
521:            grad(J23-1) = grad(J23-1) + f * dy551:            grad(J2*3-1) = grad(J2*3-1) + f * dy
522:            grad(J23)   = grad(J23)   + f * dz552:            grad(J2*3)   = grad(J2*3)   + f * dz
523:        ENDDO 553:        ENDDO 
524: !$OMP END DO554: !$OMP END DO
525: !$OMP END PARALLEL555: !$OMP END PARALLEL
526: 556: 
527:       END557:       END
528: 558: 
529: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END OF SBMBONDS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^559: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END OF SBMBONDS^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
530: 560: 
531: 561: 
532: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<562: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
533: !* SBMANGL  computes the Force due to the bond angles                   *563: !* SBMANGL  computes the Force due to the bond angles                   *
534: !* This code is modeled after how AMBER performs angle forces           *564: !* This code is modeled after how AMBER performs angle forces           *
535: !***********************************************************************565: !***********************************************************************
536: 566: 
537:       SUBROUTINE SBMANGL(grad, energy, NATOMS)567:       SUBROUTINE SBMANGL(x,y,z,grad, energy, NATOMS,IT, JT, KT, 
 568:      Q ANTC, Tk,NTA)
538:       USE KEY569:       USE KEY
539:       USE SBMDATA570:       implicit NONE
540:       IMPLICIT NONE 
541:       integer NATOMS,C1571:       integer NATOMS,C1
542:       DOUBLE PRECISION grad(3*NATOMS),energy572:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), grad(3*NATOMS),
 573:      Q energy
 574:       integer NTA
543:       LOGICAL SKIP,NOCRST575:       LOGICAL SKIP,NOCRST
544: 576: 
545:       DOUBLE PRECISION RIJ,RKJ,RIK,ANT,XIJ,YIJ,577:       DOUBLE PRECISION RIJ,RKJ,RIK,ANT,XIJ,YIJ,
546:      + ZIJ,XKJ,YKJ,ZKJ, DF578:      + ZIJ,XKJ,YKJ,ZKJ, DF
547:       DOUBLE PRECISION CT0, CT1, CT2, DA, ST, 579:       DOUBLE PRECISION CT0, CT1, CT2, DA, ST, 
548:      + CIK, CII, CKK, DT1, DT2, DT3, DT4, DT5, DT6, DT7, DT8, DT9,  STH580:      + CIK, CII, CKK, DT1, DT2, DT3, DT4, DT5, DT6, DT7, DT8, DT9,  STH
549:       DOUBLE PRECISION ONE,NEGONE581:       DOUBLE PRECISION ONE,NEGONE
550: 582: 
551:         INTEGER II583:         DOUBLE PRECISION ANTC(NTA), Tk(NTA)
552:         INTEGER I3,J3,K3,I33,J33,K33584:         INTEGER II, IT(NTA), JT(NTA), KT(NTA)
 585:         INTEGER I3, J3, K3
553:         ONE=1.0586:         ONE=1.0
554:         NEGONE=-1.0587:         NEGONE=-1.0
555: !$OMP PARALLEL PRIVATE(I3,J3,K3,I33,J33,K33,RIJ,RKJ,RIK,ANT,XIJ,YIJ,ZIJ,XKJ,YKJ,588: !$OMP PARALLEL PRIVATE(RIJ,RKJ,RIK,ANT,XIJ,YIJ,ZIJ,XKJ,YKJ,
556: !$OMP&  ZKJ,DF,CT0,CT1,CT2,DA,ST,CIK,CII,CKK,DT1,DT2,DT3,DT4,589: !$OMP&  ZKJ,DF,CT0,CT1,CT2,DA,ST,CIK,CII,CKK,DT1,DT2,DT3,DT4,
557: !$OMP&  DT5,DT6,DT7,DT8,DT9,STH) REDUCTION(+:ENERGY,grad)590: !$OMP&  DT5,DT6,DT7,DT8,DT9,STH) REDUCTION(+:ENERGY,grad)
558: !$OMP DO591: !$OMP DO
559:           DO II = 1, nTA592:           DO II = 1, nTA
560:             I3 = IT(II)593:             I3 = IT(II)
561:             J3 = JT(II)594:             J3 = JT(II)
562:             K3 = KT(II)595:             K3 = KT(II)
563:             I33=I3*3596:             XIJ = X(I3)-X(J3)
564:             J33=J3*3597:             YIJ = Y(I3)-Y(J3)
565:             K33=K3*3598:             ZIJ = Z(I3)-Z(J3)
566:             XIJ = XSBM(I3)-XSBM(J3)599:             XKJ = X(K3)-X(J3)
567:             YIJ = YSBM(I3)-YSBM(J3)600:             YKJ = Y(K3)-Y(J3)
568:             ZIJ = ZSBM(I3)-ZSBM(J3)601:             ZKJ = Z(K3)-Z(J3)
569:             XKJ = XSBM(K3)-XSBM(J3) 
570:             YKJ = YSBM(K3)-YSBM(J3) 
571:             ZKJ = ZSBM(K3)-ZSBM(J3) 
572:             RIJ = XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ602:             RIJ = XIJ*XIJ+YIJ*YIJ+ZIJ*ZIJ
573:             RKJ = XKJ*XKJ+YKJ*YKJ+ZKJ*ZKJ603:             RKJ = XKJ*XKJ+YKJ*YKJ+ZKJ*ZKJ
574:             RIK = SQRT(RIJ*RKJ)604:             RIK = SQRT(RIJ*RKJ)
575:             CT0 = (XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RIK605:             CT0 = (XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RIK
576:             CT1 = MAX(NEGONE,CT0)606:             CT1 = MAX(NEGONE,CT0)
577:             CT2 = MIN(ONE,CT1)607:             CT2 = MIN(ONE,CT1)
578:             ANT = ACOS(CT2)608:             ANT = ACOS(CT2)
579:             DA = ANT - ANTC(II)609:             DA = ANT - ANTC(II)
580:             DF = TK(II)*DA610:             DF = TK(II)*DA
581:             ST = -(DF)/SIN(ANT)611:             ST = -(DF)/SIN(ANT)
585:             CKK = STH/RKJ615:             CKK = STH/RKJ
586:             DT1 = CIK*XKJ-CII*XIJ616:             DT1 = CIK*XKJ-CII*XIJ
587:             DT2 = CIK*YKJ-CII*YIJ617:             DT2 = CIK*YKJ-CII*YIJ
588:             DT3 = CIK*ZKJ-CII*ZIJ618:             DT3 = CIK*ZKJ-CII*ZIJ
589:             DT7 = CIK*XIJ-CKK*XKJ619:             DT7 = CIK*XIJ-CKK*XKJ
590:             DT8 = CIK*YIJ-CKK*YKJ620:             DT8 = CIK*YIJ-CKK*YKJ
591:             DT9 = CIK*ZIJ-CKK*ZKJ621:             DT9 = CIK*ZIJ-CKK*ZKJ
592:             DT4 = -DT1-DT7622:             DT4 = -DT1-DT7
593:             DT5 = -DT2-DT8623:             DT5 = -DT2-DT8
594:             DT6 = -DT3-DT9624:             DT6 = -DT3-DT9
595:             grad(I33-2) = grad(I33-2)+ DT1625:             grad(I3*3-2) = grad(I3*3-2)+ DT1
596:             grad(I33-1) = grad(I33-1)+ DT2626:             grad(I3*3-1) = grad(I3*3-1)+ DT2
597:             grad(I33)   = grad(I33)  + DT3627:             grad(I3*3)   = grad(I3*3)  + DT3
598:             grad(J33-2) = grad(J33-2)+ DT4628:             grad(J3*3-2) = grad(J3*3-2)+ DT4
599:             grad(J33-1) = grad(J33-1)+ DT5629:             grad(J3*3-1) = grad(J3*3-1)+ DT5
600:             grad(J33)   = grad(J33)  + DT6630:             grad(J3*3)   = grad(J3*3)  + DT6
601:             grad(K33-2) = grad(K33-2)+ DT7631:             grad(K3*3-2) = grad(K3*3-2)+ DT7
602:             grad(K33-1) = grad(K33-1)+ DT8632:             grad(K3*3-1) = grad(K3*3-1)+ DT8
603:             grad(K33)   = grad(K33)  + DT9633:             grad(K3*3)   = grad(K3*3)  + DT9
604:             ENERGY = ENERGY + TK(II)*(ANTC(II)- ANT)**2/2.0634:             ENERGY = ENERGY + TK(II)*(ANTC(II)- ANT)**2/2.0
605:           END DO635:           END DO
606: !$OMP END DO636: !$OMP END DO
607: !$OMP END PARALLEL637: !$OMP END PARALLEL
608: 638: 
609:        RETURN639:        RETURN
610:        END640:        END
611: 641: 
612: !^^^^^^^^^^^^^^^^^^^^^^^^End of SBMANGL^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^642: !^^^^^^^^^^^^^^^^^^^^^^^^End of SBMANGL^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
613: 643: 
614: 644: 
615: 645: 
616: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<646: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
617: !* SBMdihedral computes the dihedral angles and the forces due to them *647: !* SBMdihedral computes the dihedral angles and the forces due to them *
618: !**********************************************************************648: !**********************************************************************
619: 649: 
620:       SUBROUTINE SBMdihedral(grad,energy, NATOMS)650:       SUBROUTINE SBMdihedral(x,y,z,grad, energy, NATOMS,IP,JP,KP,LP,PK,
 651:      Q PHITYPE,PHISBM,NPA)
621:       USE KEY652:       USE KEY
622:       USE SBMDATA 
623:       implicit NONE653:       implicit NONE
624:       integer I, N, J, NATOMS, II654:       integer I, N, J, NATOMS, NPA, II
625:       DOUBLE PRECISION grad(3*NATOMS),energy655:       DOUBLE PRECISION x(NATOMS),y(NATOMS),z(NATOMS),
 656:      Q grad(3*NATOMS),energy
 657: 
 658: 
 659:       DOUBLE PRECISION PK(NPA),PHISBM(NPA)
 660:       INTEGER IP(NPA), JP(NPA), KP(NPA), LP(NPA) 
 661:       INTEGER PHITYPE(NPA)
626: 662: 
627:       double precision lfac663:       double precision lfac
628:       integer I3,J3,K3,L3,C1,I33,J33,K33,L33664:       integer I3, J3, K3, L3,C1
629:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,665:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,
630:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,666:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,
631:      + DZ, GX,GY,GZ,CT,CPHI,667:      + DZ, GX,GY,GZ,CT,CPHI,
632:      + SPHI,Z1, Z2,FXI,FYI,FZI,668:      + SPHI,Z1, Z2,FXI,FYI,FZI,
633:      + FXJ,FYJ,FZJ, FXK,FYK,FZK,669:      + FXJ,FYJ,FZJ, FXK,FYK,FZK,
634:      + FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,ftem,CT0,CT1,AP0,AP1,670:      + FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,ftem,CT0,CT1,AP0,AP1,
635:      + Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,671:      + Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,
636:      + S,HGoverG,FGoverG,A1,A3672:      + S,HGoverG,FGoverG,A1,A3
637: 673: 
 674: 
 675: C
638:       double precision  TM24,TM06,tenm3,zero,one,NEGONE,two,four,six,twelve676:       double precision  TM24,TM06,tenm3,zero,one,NEGONE,two,four,six,twelve
639: 677: 
640:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,678:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,
641:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z679:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z
642: 680: 
643:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/681:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/
644: 682: 
645:       double precision pi683:       double precision pi
646:       pi = 3.14159265358979323846264338327950288419716939937510684:       pi = 3.14159265358979323846264338327950288419716939937510
647:       ONE=1.0685:       ONE=1.0
648:       NEGONE=-1.0686:       NEGONE=-1.0
649: 687: 
650: !$OMP PARALLEL PRIVATE(I,I3,J3,K3,L3,I33,J33,K33,L33,II,XIJ,YIJ,ZIJ,XKJ,YKJ, 688: !$OMP PARALLEL PRIVATE(I,I3,J3,K3,L3,II,XIJ,YIJ,ZIJ,XKJ,YKJ, 
651: !$OMP& ZKJ,XKL,YKL,ZKL,RIJ,RKJ,RKL,DX,DY, 689: !$OMP& ZKJ,XKL,YKL,ZKL,RIJ,RKJ,RKL,DX,DY, 
652: !$OMP& DZ, GX,GY,GZ,CT,CPHI,690: !$OMP& DZ, GX,GY,GZ,CT,CPHI,
653: !$OMP& SPHI,Z1, Z2,FXI,FYI,FZI,691: !$OMP& SPHI,Z1, Z2,FXI,FYI,FZI,
654: !$OMP& FXJ,FYJ,FZJ, FXK,FYK,FZK,692: !$OMP& FXJ,FYJ,FZJ, FXK,FYK,FZK,
655: !$OMP& FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,CT0,CT1,AP0,AP1,693: !$OMP& FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,CT0,CT1,AP0,AP1,
656: !$OMP& Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,694: !$OMP& Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,
657: !$OMP& S,HGoverG,FGoverG,A1,A3,TT1,TT2,TT3,TT4,TT1X,TT1Y,TT1Z,TT2X,695: !$OMP& S,HGoverG,FGoverG,A1,A3,TT1,TT2,TT3,TT4,TT1X,TT1Y,TT1Z,TT2X,
658: !$OMP& TT2Y,TT2Z,TT3X,TT3Y,TT3Z,TT4X,TT4Y,TT4Z) REDUCTION(+:energy,grad)696: !$OMP& TT2Y,TT2Z,TT3X,TT3Y,TT3Z,TT4X,TT4Y,TT4Z) REDUCTION(+:energy,grad)
659: !!697: !!
660: !$OMP DO698: !$OMP DO
661:           DO II = 1,nPA699:           DO II = 1,nPA
662: 700: 
663:             I3 = IP(II)701:             I3 = IP(II)
664:             J3 = JP(II)702:             J3 = JP(II)
665:             K3 = KP(II)703:             K3 = KP(II)
666:             L3 = LP(II)704:             L3 = LP(II)
667:             I33=I3*3705:             XIJ = X(I3)-X(J3)
668:             J33=J3*3706:             YIJ = Y(I3)-Y(J3)
669:             K33=K3*3707:             ZIJ = Z(I3)-Z(J3)
670:             L33=L3*3708:             XKJ = X(K3)-X(J3)
671:             XIJ = XSBM(I3)-XSBM(J3)709:             YKJ = Y(K3)-Y(J3)
672:             YIJ = YSBM(I3)-YSBM(J3)710:             ZKJ = Z(K3)-Z(J3)
673:             ZIJ = ZSBM(I3)-ZSBM(J3) 
674:             XKJ = XSBM(K3)-XSBM(J3) 
675:             YKJ = YSBM(K3)-YSBM(J3) 
676:             ZKJ = ZSBM(K3)-ZSBM(J3) 
677:             RKJ = sqrt(XKJ**2+YKJ**2+ZKJ**2)711:             RKJ = sqrt(XKJ**2+YKJ**2+ZKJ**2)
678:             XKL = XSBM(K3)-XSBM(L3)712:             XKL = X(K3)-X(L3)
679:             YKL = YSBM(K3)-YSBM(L3)713:             YKL = Y(K3)-Y(L3)
680:             ZKL = ZSBM(K3)-ZSBM(L3)                                  714:             ZKL = Z(K3)-Z(L3)                                  
681:             FGoverG=-(XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RKJ715:             FGoverG=-(XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RKJ
682:             HGoverG=(XKL*XKJ+YKL*YKJ+ ZKL*ZKJ)/RKJ716:             HGoverG=(XKL*XKJ+YKL*YKJ+ ZKL*ZKJ)/RKJ
683: C DX is the M vector and G is the N vector717: C DX is the M vector and G is the N vector
684:             DX = YIJ*ZKJ-ZIJ*YKJ718:             DX = YIJ*ZKJ-ZIJ*YKJ
685:             DY = ZIJ*XKJ-XIJ*ZKJ719:             DY = ZIJ*XKJ-XIJ*ZKJ
686:             DZ = XIJ*YKJ-YIJ*XKJ720:             DZ = XIJ*YKJ-YIJ*XKJ
687:             GX = ZKJ*YKL-YKJ*ZKL721:             GX = ZKJ*YKL-YKJ*ZKL
688:             GY = XKJ*ZKL-ZKJ*XKL722:             GY = XKJ*ZKL-ZKJ*XKL
689:             GZ = YKJ*XKL-XKJ*YKL723:             GZ = YKJ*XKL-XKJ*YKL
690:             FXI = SQRT(DX*DX+DY*DY+DZ*DZ)724:             FXI = SQRT(DX*DX+DY*DY+DZ*DZ)
746:           TT1Z=TT1*DZ780:           TT1Z=TT1*DZ
747:           TT2X=TT2*DX781:           TT2X=TT2*DX
748:           TT2Y=TT2*DY782:           TT2Y=TT2*DY
749:           TT2Z=TT2*DZ783:           TT2Z=TT2*DZ
750:           TT3X=TT3*GX784:           TT3X=TT3*GX
751:           TT3Y=TT3*GY785:           TT3Y=TT3*GY
752:           TT3Z=TT3*GZ786:           TT3Z=TT3*GZ
753:           TT4X=TT4*GX787:           TT4X=TT4*GX
754:           TT4Y=TT4*GY788:           TT4Y=TT4*GY
755:           TT4Z=TT4*GZ789:           TT4Z=TT4*GZ
756:           grad(I33-2) =  grad(I33-2)  + TT1X  790:           grad(I3*3-2) =  grad(I3*3-2)  + TT1X  
757:           grad(I33-1) =  grad(I33-1)  + TT1Y791:           grad(I3*3-1) =  grad(I3*3-1)  + TT1Y
758:           grad(I33)   =  grad(I33)    + TT1Z792:           grad(I3*3)   =  grad(I3*3)    + TT1Z
759:           grad(J33-2) =  grad(J33-2)  - TT1X - TT2X - TT3X793:           grad(J3*3-2) =  grad(J3*3-2)  - TT1X - TT2X - TT3X
760:           grad(J33-1) =  grad(J33-1)  - TT1Y - TT2Y - TT3Y794:           grad(J3*3-1) =  grad(J3*3-1)  - TT1Y - TT2Y - TT3Y
761:           grad(J33)   =  grad(J33)    - TT1Z - TT2Z - TT3Z795:           grad(J3*3)   =  grad(J3*3)    - TT1Z - TT2Z - TT3Z
762:           grad(K33-2) =  grad(K33-2)  + TT2X + TT3X - TT4X796:           grad(K3*3-2) =  grad(K3*3-2)  + TT2X + TT3X - TT4X
763:           grad(K33-1) =  grad(K33-1)  + TT2Y + TT3Y - TT4Y797:           grad(K3*3-1) =  grad(K3*3-1)  + TT2Y + TT3Y - TT4Y
764:           grad(K33)   =  grad(K33)    + TT2Z + TT3Z - TT4Z798:           grad(K3*3)   =  grad(K3*3)    + TT2Z + TT3Z - TT4Z
765:           grad(L33-2) =  grad(L33-2)  + TT4X799:           grad(L3*3-2) =  grad(L3*3-2)  + TT4X
766:           grad(L33-1) =  grad(L33-1)  + TT4Y800:           grad(L3*3-1) =  grad(L3*3-1)  + TT4Y
767:           grad(L33)   =  grad(L33)    + TT4Z801:           grad(L3*3)   =  grad(L3*3)    + TT4Z
768:         END DO802:         END DO
769: !$OMP END DO 803: !$OMP END DO 
770: !$OMP END PARALLEL804: !$OMP END PARALLEL
771: 805: 
772:           END806:           END
773: 807: 
774: 808: 
775: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^809: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^
776: 810: 
777: 811: 
778: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<812: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
779: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *813: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *
780: !* 10-12 or 6-12 potential                                              *814: !* 10-12 or 6-12 potential                                              *
781: !***********************************************************************815: !***********************************************************************
782: 816: 
783:       subroutine SBMcontacts(grad,energy,NATOMS)817:       subroutine SBMcontacts(x,y,z,grad,energy,
 818:      Q NATOMS,IC,JC,CONCOEF,NC,CONTACTTYPE)
784:       USE KEY819:       USE KEY
785:       USE SBMDATA 
786:       implicit NONE820:       implicit NONE
787:       integer I,J,NATOMS821:       integer I, NATOMS,NC
788: 822: 
789:       DOUBLE PRECISION grad(3*NATOMS), energy823:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
 824:      Q grad(3*NATOMS), energy
790:       DOUBLE PRECISION dx,dy,dz,CO1,CO2,CO3,CO4,CO5,CO6,E1,G1,G2,E1D,G1D,G2D,ET825:       DOUBLE PRECISION dx,dy,dz,CO1,CO2,CO3,CO4,CO5,CO6,E1,G1,G2,E1D,G1D,G2D,ET
791: 826: 
792:       integer C1,C2,C13,C23 827:       integer C1, C2 
793:       DOUBLE PRECISION  r2,rm2,rm10,rm12,rm6,f_over_r,r1 828:       DOUBLE PRECISION  r2,rm2,rm10,rm12,rm6,f_over_r,r1 
794: 829: 
795: ! CONTACTTYPE830:         DOUBLE PRECISION CONCOEF(NC*6)
 831:         INTEGER IC(NC), JC(NC),CONTACTTYPE(NC)
 832: 
796: ! type 1 is 6-12833: ! type 1 is 6-12
797: ! type 2 is 12-12834: ! type 2 is 12-12
798: ! type 5 as gaussian, no ex-vol835: ! implement type 3 as gaussian
799: ! type 6 as gaussian, w ex-vol836: ! implement type 4 as dual-basin gaussian
800: ! type 7 as dual gaussian837: ! question: Should type be part of the interaction, or should they be organized
801: !$OMP PARALLEL PRIVATE(I,CO1,CO2,CO3,CO4,CO5,CO6,C1,C2,C13,C23,DX,DY,DZ,838: ! into separate lists?
 839: ! type 1 is 6-12 interaction
 840: !$OMP PARALLEL PRIVATE(I,CO1,CO2,CO3,CO4,CO5,CO6,C1,C2,DX,DY,DZ,
802: !$OMP& R2,RM2,RM6,RM10,RM12,R1,E1,G1,G2,E1D,G1D,G2D,F_OVER_R,ET) REDUCTION(+:ENERGY,grad)841: !$OMP& R2,RM2,RM6,RM10,RM12,R1,E1,G1,G2,E1D,G1D,G2D,F_OVER_R,ET) REDUCTION(+:ENERGY,grad)
803: !$OMP DO842: !$OMP DO
804:           do i=1, NC843:           do i=1, NC
805:             C1 = IC(i)844:             C1 = IC(i)
806:             C2 = JC(i)845:             C2 = JC(i)
807:             dx = XSBM(C1) - XSBM(C2)846:             dx = X(C1) - X(C2)
808:             dy = YSBM(C1) - YSBM(C2)847:             dy = Y(C1) - Y(C2)
809:             dz = ZSBM(C1) - ZSBM(C2)848:             dz = Z(C1) - Z(C2)
810:             r2 = dx**2 + dy**2 + dz**2849:             r2 = dx**2 + dy**2 + dz**2
811:             J=I*6 
812:             C13=C1*3 
813:             C23=C2*3 
814:             if(CONTACTTYPE(I) .eq. 1)then850:             if(CONTACTTYPE(I) .eq. 1)then
815:               ! type 1 is 6-12 interaction851:               ! type 1 is 6-12 interaction
816:               rm2 = 1.0/r2852:               rm2 = 1.0/r2
817:               rm6 = rm2**3853:               rm6 = rm2**3
818:               CO1=CONCOEF(J-5)854:               CO1=CONCOEF(I*6-5)
819:               CO2=CONCOEF(J-4)855:               CO2=CONCOEF(I*6-4)
820:               ENERGY = ENERGY + rm6*(CO1*rm6-CO2)856:               ENERGY = ENERGY + rm6*(CO1*rm6-CO2)
821:               f_over_r = -12.0*rm6*(CO1*rm6-CO2/2.0)*rm2857:               f_over_r = -12.0*rm6*(CO1*rm6-CO2/2.0)*rm2
822:             elseif(CONTACTTYPE(I) .eq. 2)then858:             elseif(CONTACTTYPE(I) .eq. 2)then
823:               ! type 2 is 10-12 interaction859:               ! type 2 is 10-12 interaction
824:               rm2 = 1.0/r2860:               rm2 = 1.0/r2
825:               rm10 = rm2**5861:               rm10 = rm2**5
826:               CO1=CONCOEF(J-5)862:               CO1=CONCOEF(I*6-5)
827:               CO2=CONCOEF(J-4)863:               CO2=CONCOEF(I*6-4)
828:               ENERGY = ENERGY + rm10*(CO1*rm2-CO2)864:               ENERGY = ENERGY + rm10*(CO1*rm2-CO2)
829:               f_over_r = -60.0*rm10*(CO1*rm2/5.0-CO2/6.0)*rm2865:               f_over_r = -60.0*rm10*(CO1*rm2/5.0-CO2/6.0)*rm2
830:             elseif(CONTACTTYPE(I) .eq. 5)then866:             elseif(CONTACTTYPE(I) .eq. 5)then
831:               ! type 5 is a simple gaussian interaction (no excluded867:               ! type 5 is a simple gaussian interaction (no excluded
832:               ! volume)868:               ! volume)
833:               r1 = sqrt(r2)869:               r1 = sqrt(r2)
834:               CO1=CONCOEF(J-5) !A870:               CO1=CONCOEF(I*6-5) !A
835:               CO2=CONCOEF(J-4) !mu871:               CO2=CONCOEF(I*6-4) !mu
836:               CO3=CONCOEF(J-3) !1/(2*sigma**2)872:               CO3=CONCOEF(I*6-3) !1/(2*sigma**2)
837:               E1= -CO1*exp(-CO3*(r1-CO2)**2)873:               E1= -CO1*exp(-CO3*(r1-CO2)**2)
838:               ENERGY = ENERGY + E1874:               ENERGY = ENERGY + E1
839:               f_over_r = -E1*2*CO3*(r1-CO2)/r1875:               f_over_r = -E1*2*CO3*(r1-CO2)/r1
840:             elseif(CONTACTTYPE(I) .eq. 6)then876:             elseif(CONTACTTYPE(I) .eq. 6)then
841:               ! type 6 is a gaussian interaction with excluded877:               ! type 6 is a gaussian interaction with excluded
842:               ! volume included878:               ! volume included
843:               r1 = sqrt(r2)879:               r1 = sqrt(r2)
844:               rm2 = 1.0/r2880:               rm2 = 1.0/r2
845:               rm12 = rm2**6881:               rm12 = rm2**6
846:               CO1=CONCOEF(J-5) !A882:               CO1=CONCOEF(I*6-5) !A
847:               CO2=CONCOEF(J-4) !mu883:               CO2=CONCOEF(I*6-4) !mu
848:               CO3=CONCOEF(J-3) !1/(2*sigma**2)884:               CO3=CONCOEF(I*6-3) !1/(2*sigma**2)
849:               CO4=CONCOEF(J-2) !a/A885:               CO4=CONCOEF(I*6-2) !a/A
850:               E1 =CO1*(1+CO4*rm12)886:               E1 =CO1*(1+CO4*rm12)
851:               G1 =1-exp(-CO3*(r1-CO2)**2)887:               G1 =1-exp(-CO3*(r1-CO2)**2)
852:               E1D=-12*CO1*CO4*rm12*rm2888:               E1D=-12*CO1*CO4*rm12*rm2
853:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)889:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)
854:               ENERGY = ENERGY + (E1*G1-CO1)890:               ENERGY = ENERGY + (E1*G1-CO1)
855:               f_over_r = E1*G1D + G1*E1D 891:               f_over_r = E1*G1D + G1*E1D 
856:             elseif(CONTACTTYPE(I) .eq. 7)then892:             elseif(CONTACTTYPE(I) .eq. 7)then
857:               ! type 7 is a double gaussian interaction with excluded893:               ! type 7 is a double gaussian interaction with excluded
858:               ! volume included894:               ! volume included
859:               rm2 = 1.0/r2895:               rm2 = 1.0/r2
860:               rm12 = rm2**6896:               rm12 = rm2**6
861:               r1 = sqrt(r2)897:               r1 = sqrt(r2)
862:               CO1=CONCOEF(J-5) !A898:               CO1=CONCOEF(I*6-5) !A
863:               CO2=CONCOEF(J-4) !mu1899:               CO2=CONCOEF(I*6-4) !mu1
864:               CO3=CONCOEF(J-3) !1/(2*sigma1**2)900:               CO3=CONCOEF(I*6-3) !1/(2*sigma1**2)
865:               CO4=CONCOEF(J-2) !mu2901:               CO4=CONCOEF(I*6-2) !mu2
866:               CO5=CONCOEF(J-1) !1/(2*sigma2**2)902:               CO5=CONCOEF(I*6-1) !1/(2*sigma2**2)
867:               CO6=CONCOEF(J)   ! a/A (excluded volume amplitude/gaussian amplitude)903:               CO6=CONCOEF(I*6)   ! a/A (excluded volume amplitude/gaussian amplitude)
868:               E1=1+CO6*rm12904:               E1=1+CO6*rm12
869:               G1=1-exp(-CO3*(r1-CO2)**2)905:               G1=1-exp(-CO3*(r1-CO2)**2)
870:               G2=1-exp(-CO5*(r1-CO4)**2)906:               G2=1-exp(-CO5*(r1-CO4)**2)
871:               E1D= -12*CO6*rm12*rm2907:               E1D= -12*CO6*rm12*rm2
872:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)908:               G1D=CO3*2*(1-CO2/r1)*exp(-CO3*(r1-CO2)**2)
873:               G2D=CO5*2*(1-CO4/r1)*exp(-CO5*(r1-CO4)**2)909:               G2D=CO5*2*(1-CO4/r1)*exp(-CO5*(r1-CO4)**2)
874:               ENERGY = ENERGY + CO1*( E1 * G1 * G2-1) 910:               ENERGY = ENERGY + CO1*( E1 * G1 * G2-1) 
875:               f_over_r = CO1*(E1D*G1*G2 911:               f_over_r = CO1*(E1D*G1*G2 
876:      Q        +  E1 * G1D * G2912:      Q        +  E1 * G1D * G2
877:      Q        +  E1 * G1 * G2D)913:      Q        +  E1 * G1 * G2D)
878:             else914:             else
879:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)915:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)
880:               STOP916:               STOP
881:             endif917:             endif
882: 918: 
883:             grad(C13-2) = grad(C13-2) + f_over_r * dx919:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
884:             grad(C13-1) = grad(C13-1) + f_over_r * dy920:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
885:             grad(C13)   = grad(C13)   + f_over_r * dz921:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz
886:             grad(C23-2) = grad(C23-2) - f_over_r * dx922:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
887:             grad(C23-1) = grad(C23-1) - f_over_r * dy923:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
888:             grad(C23)   = grad(C23)   - f_over_r * dz924:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz
889:           enddo925:           enddo
890: !$OMP END DO926: !$OMP END DO
891: !$OMP END PARALLEL927: !$OMP END PARALLEL
892: 928: 
893: 929: 
894:       end930:       end
895: 931: 
896: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^932: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
897: 933: 
898: 934: 
899: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<935: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
900: !* SBMNonContacts computes the forces due to non native contacts      *936: !* SBMNonContacts computes the forces due to non native contacts      *
901: !**********************************************************************937: !**********************************************************************
902: 938: 
903:         subroutine SBMnoncontacts(grad,energy,NATOMS)939:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,ATOMTYPES,
 940:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,
 941:      Q NCswitch,NCcut,STT,SA,SB,SC)
904:         USE KEY942:         USE KEY
905:         USE SBMDATA 
906:         implicit NONE943:         implicit NONE
907:         integer I,I3,N,J,AN,NATOMS944:         integer I,N,J,AN,NATOMS,ATOMTYPES(NATOMS),MAXATOMTYPES
908: 945: 
909:         DOUBLE PRECISION grad(3*NATOMS), energy946:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
910:        DOUBLE PRECISION SAT,SBT,SCT,STTT947:      Q  grad(3*NATOMS), energy
911: 948:        DOUBLE PRECISION STT,SA,SB,SC,SAT,SBT,SCT,STTT
912:         integer C1,C2,C13,C23,DINDEX,ii,jj,kk,k,l,iii,jjj,POS949:        DIMENSION STT(MAXATOMTYPES*MAXATOMTYPES)
913:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r 950:        DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
914:         INTEGER NTHREADS951:        DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
 952:        DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
 953: 
 954: 
 955:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj,POS
 956:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
 957:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,MAXSEPSYS
 958:         INTEGER NNCINC
 959:         DIMENSION NNCINC(NATOMS*MAXSEP)
915: !$    INTEGER  OMP_GET_NUM_THREADS,960: !$    INTEGER  OMP_GET_NUM_THREADS,
916: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS961: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
917: 962: 
 963:         INTEGER NNEXL1(NEXCLUSIONS)
 964:         INTEGER NNEXL2(NEXCLUSIONS)
918:         DOUBLE PRECISION dx,dy,dz965:         DOUBLE PRECISION dx,dy,dz
919:         integer tempN, alpha966:         integer tempN, alpha
920:         double precision Rdiff,Vfunc,Ffunc967:         double precision Rdiff,Vfunc,Ffunc
921:         double precision Rcut2,Rswitch2968:         double precision Rcut2,Rswitch2
922:         ! Ngrid is the number of atoms in that grid point969:         ! Ngrid is the number of atoms in that grid point
923:         ! grid is the array of atoms in each grid970:         ! grid is the array of atoms in each grid
924:         integer Ngrid,grid,maxgrid,maxpergrid971:         integer Ngrid,grid,maxgrid,maxpergrid
925:         ! number of atoms per grid, max972:         ! number of atoms per grid, max
926:         parameter (maxpergrid=50)973:         parameter (maxpergrid=50)
927:         ! dimensions of grid974:         ! dimensions of grid
942:         minX=10000000989:         minX=10000000
943:         minY=10000000990:         minY=10000000
944:         minZ=10000000991:         minZ=10000000
945:         992:         
946:         maxX=-10000000993:         maxX=-10000000
947:         maxY=-10000000994:         maxY=-10000000
948:         maxZ=-10000000995:         maxZ=-10000000
949: 996: 
950:         do i=1,NATOMS997:         do i=1,NATOMS
951: 998: 
952:            if(XSBM(i) .lt. minX)then999:            if(X(i) .gt. maxX)then
953:                 minX=XSBM(i)1000:                 maxX=X(i)
954:            elseif(XSBM(i) .gt. maxX)then1001:            elseif(Y(i) .gt. maxY)then
955:                 maxX=XSBM(i)1002:                 maxY=Y(i)
956:            endif 
957:            if(YSBM(i) .lt. minY)then 
958:                 minY=YSBM(i) 
959:            elseif(YSBM(i) .gt. maxY)then 
960:                 maxY=YSBM(i) 
961:            endif1003:            endif
962:            if(ZSBM(i) .lt. minZ)then1004: 
963:                 minZ=ZSBM(i)1005:            if(Z(i) .gt. maxZ)then
964:            elseif(ZSBM(i) .gt. maxZ)then1006:                 maxZ=Z(i)
965:                 maxZ=ZSBM(i)1007:            elseif(X(i) .lt. minX)then
 1008:                 minX=X(i)
966:            endif1009:            endif
967: 1010: 
 1011:            if(Y(i) .lt. minY)then
 1012:                 minY=Y(i)
 1013:            elseif(Z(i) .lt. minZ)then
 1014:                 minZ=Z(i)
 1015:            endif
968:         enddo1016:         enddo
969: 1017: 
970:         maxgridX=int((maxX-minX)/gridsize)+11018:         maxgridX=int((maxX-minX)/gridsize)+1
971:         maxgridY=int((maxY-minY)/gridsize)+11019:         maxgridY=int((maxY-minY)/gridsize)+1
972:         maxgridZ=int((maxZ-minZ)/gridsize)+11020:         maxgridZ=int((maxZ-minZ)/gridsize)+1
973: 1021: 
974:         if(maxgridX .ge. maxgrid .or. 1022:         if(maxgridX .ge. maxgrid .or. 
975:      Q  maxgridY .ge. maxgrid .or.1023:      Q  maxgridY .ge. maxgrid .or.
976:      Q  maxgridZ .ge. maxgrid )then1024:      Q  maxgridZ .ge. maxgrid )then
977:         write(*,*) 'ERROR: system got too big for grid searching...'1025:         write(*,*) 'ERROR: system got too big for grid searching...'
982: !!  Add a second grid that only includes atoms that are charged1030: !!  Add a second grid that only includes atoms that are charged
983: 1031: 
984:         do i=1,maxgrid1032:         do i=1,maxgrid
985:          do j=1,maxgrid1033:          do j=1,maxgrid
986:           do k=1,maxgrid1034:           do k=1,maxgrid
987:                 Ngrid(i,j,k)=01035:                 Ngrid(i,j,k)=0
988:           enddo1036:           enddo
989:          enddo1037:          enddo
990:         enddo1038:         enddo
991:         do i=1,NATOMS1039:         do i=1,NATOMS
992:                 Xgrid=int((XSBM(i)-minX)/gridsize)+11040: 
993:                 Ygrid=int((YSBM(i)-minY)/gridsize)+11041:                 Xgrid=int((X(i)-minX)/gridsize)+1
994:                 Zgrid=int((ZSBM(i)-minZ)/gridsize)+11042:                 Ygrid=int((Y(i)-minY)/gridsize)+1
 1043:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
995:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+11044:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
996:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then1045:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
997:                         write(*,*) 'ERROR: too many atoms in a grid'1046:                         write(*,*) 'ERROR: too many atoms in a grid'
998:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid1047:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
999:                         STOP1048:                         STOP
1000:                 endif1049:                 endif
1001:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i1050:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1002:         enddo1051:         enddo
1003: 1052: 
1004: 1053: 
1005:            tempN=01054:            tempN=0
1006: 1055: 
1007: ! add a second loop that goes over charged atoms1056: ! add a second loop that goes over charged atoms
1008: 1057: 
1009: !$OMP PARALLEL1058: !$OMP PARALLEL
1010: !$OMP& PRIVATE(i,I3,ii,jj,kk,jjj,C1,C2,C13,C23,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,1059: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,
1011: !$OMP& Xgrid,Ygrid,Zgrid,TID,DINDEX,POS,SAT,SBT,SCT,STTT)  1060: !$OMP& Xgrid,Ygrid,Zgrid,TID,DINDEX,POS,SAT,SBT,SCT,STTT)  
1012: !$OMP& REDUCTION(+:ENERGY,grad)1061: !$OMP& REDUCTION(+:ENERGY,grad)
1013:         TID=01062:         TID=0
1014:         NTHREADS=1;1063:         NTHREADS=1;
1015: !$      TID = OMP_GET_THREAD_NUM()1064: !$      TID = OMP_GET_THREAD_NUM()
1016: !$      NTHREADS = OMP_GET_NUM_THREADS()1065: !$      NTHREADS = OMP_GET_NUM_THREADS()
1017:         do i=1+TID,NATOMS,NTHREADS1066:         do i=1+TID,NATOMS,NTHREADS
1018:           I3=I*31067: 
1019:           Xgrid=int((XSBM(I)-minX)/gridsize)+11068:           Xgrid=int((X(I)-minX)/gridsize)+1
1020:           Ygrid=int((YSBM(I)-minY)/gridsize)+11069:           Ygrid=int((Y(I)-minY)/gridsize)+1
1021:           Zgrid=int((ZSBM(I)-minZ)/gridsize)+11070:           Zgrid=int((Z(I)-minZ)/gridsize)+1
1022:           do ii=XGRID-1,XGRID+11071:           do ii=XGRID-1,XGRID+1
1023:            do jj=YGRID-1,YGRID+11072:            do jj=YGRID-1,YGRID+1
1024:             do kk=ZGRID-1,ZGRID+11073:             do kk=ZGRID-1,ZGRID+1
1025:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1074:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1026:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1075:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1027:               do jjj=1,Ngrid(ii,jj,kk)1076:               do jjj=1,Ngrid(ii,jj,kk)
1028:                C2=grid(ii,jj,kk,jjj)1077:                C2=grid(ii,jj,kk,jjj)
1029:                C23=C2*3 
1030:                DINDEX=C2-I1078:                DINDEX=C2-I
1031:                if(DINDEX .GT. 0)THEN1079:                if(DINDEX .GT. 0)THEN
1032:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX) .eq. 0) )then1080:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX) .eq. 0) )then
1033:                  if(SBMRBG(I) .NE. SBMRBG(C2))THEN1081:                 !if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS )  )then
1034:                   dx = XSBM(I) - XSBM(C2)1082:                  dx = X(I) - X(C2)
1035:                   dy = YSBM(I) - YSBM(C2)1083:                  dy = Y(I) - Y(C2)
1036:                   dz = ZSBM(I) - ZSBM(C2)1084:                  dz = Z(I) - Z(C2)
1037:                   r2 = dx**2 + dy**2 + dz**21085:                  r2 = dx**2 + dy**2 + dz**2
1038:                   if(r2 .le. Rcut2)then1086:                  if(r2 .le. Rcut2)then
1039:                    POS=(ATOMTYPES(I)-1)*NATOMTYPES+ATOMTYPES(C2)1087:                   POS=(ATOMTYPES(I)-1)*MAXATOMTYPES+ATOMTYPES(C2)
1040:                    STTT=STT(POS)1088:                   STTT=STT(POS)
1041:                    SCT=SC(POS)1089:                   SCT=SC(POS)
1042:                    rm2 = 1/r21090:                   rm2 = 1/r2
1043:                    rm14 = rm2**71091:                   rm14 = rm2**7
1044:                    energy=energy+STTT*rm2**6+SCT1092:                   energy=energy+STTT*rm2**6+SCT
1045:                    f_over_r=-STTT*12.0*rm141093:                   f_over_r=-STTT*12.0*rm14
1046:                    RD1=sqrt(r2)-NCswitch1094:                   RD1=sqrt(r2)-NCswitch
1047:                    if(r2 .gt. Rswitch2)then1095:                   if(r2 .gt. Rswitch2)then
1048:                     SAT=SA(POS)1096:                    SAT=SA(POS)
1049:                     SBT=SB(POS)1097:                    SBT=SB(POS)
1050:                     f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)1098:                    f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)
1051:                     energy=energy+SAT/3.0*RD1**3+SBT/4.0*RD1**41099:                    energy=energy+SAT/3.0*RD1**3+SBT/4.0*RD1**4
1052:                    elseif(r2 .lt. Rswitch2)then1100:                   elseif(r2 .lt. Rswitch2)then
1053:                   ! normal repulsive term1101:                  ! normal repulsive term
1054:                    else1102:                   else
1055:                   ! things should have fallen in one of the previous two...1103:                  ! things should have fallen in one of the previous two...
1056:                     write(*,*) 'something went wrong with switching function'1104:                    write(*,*) 'something went wrong with switching function'
1057:                    endif1105:                   endif
1058:                    grad(I3-2) = grad(I3-2) + f_over_r * dx1106:                   grad(I*3-2) = grad(I*3-2) + f_over_r * dx
1059:                    grad(I3-1) = grad(I3-1) + f_over_r * dy1107:                   grad(I*3-1) = grad(I*3-1) + f_over_r * dy
1060:                    grad(I3)   = grad(I3)   + f_over_r * dz1108:                   grad(I*3)   = grad(I*3)   + f_over_r * dz
1061:                    grad(C23-2) = grad(C23-2) - f_over_r * dx1109:                   grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
1062:                    grad(C23-1) = grad(C23-1) - f_over_r * dy1110:                   grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
1063:                    grad(C23)   = grad(C23)   - f_over_r * dz1111:                   grad(C2*3)   = grad(C2*3)   - f_over_r * dz
1064:                   ENDIF 
1065:                   endif1112:                   endif
1066:                 endif1113:                 endif
1067:                endif1114:                endif
1068:               enddo1115:               enddo
1069:              endif1116:              endif
1070:             enddo1117:             enddo
1071:            enddo1118:            enddo
1072:           enddo1119:           enddo
1073:          enddo1120:          enddo
1074: 1121: 
1075: !$OMP DO1122: !$OMP DO
1076:        DO I=1,NEXCLUSIONS1123:        DO I=1,NEXCLUSIONS
1077:          C1=NNEXL1(I) 1124:          C1=NNEXL1(I) 
1078:          C2=NNEXL2(I)1125:          C2=NNEXL2(I)
1079:          C13=C1*31126:          dx = X(C1) - X(C2)
1080:          C23=C2*31127:          dy = Y(C1) - Y(C2)
1081:          dx = XSBM(C1) - XSBM(C2)1128:          dz = Z(C1) - Z(C2)
1082:          dy = YSBM(C1) - YSBM(C2) 
1083:          dz = ZSBM(C1) - ZSBM(C2) 
1084:          r2 = dx**2 + dy**2 + dz**21129:          r2 = dx**2 + dy**2 + dz**2
1085:          if(r2 .le. Rcut2)then1130:          if(r2 .le. Rcut2)then
1086:            POS=(ATOMTYPES(C1)-1)*NATOMTYPES+ATOMTYPES(C2)1131:            POS=(ATOMTYPES(C1)-1)*MAXATOMTYPES+ATOMTYPES(C2)
1087:            STTT=STT(POS)1132:            STTT=STT(POS)
1088:            SCT=SC(POS)1133:            SCT=SC(POS)
1089:  1134:  
1090:            rm2 = 1/r21135:            rm2 = 1/r2
1091:            rm14 = rm2**71136:            rm14 = rm2**7
1092:             energy=energy-STTT*rm2**6-SCT1137:             energy=energy-STTT*rm2**6-SCT
1093:             f_over_r=-STTT*12.0*rm141138:             f_over_r=-STTT*12.0*rm14
1094:             RD1=sqrt(r2)-NCswitch1139:             RD1=sqrt(r2)-NCswitch
1095:           if(r2 .gt. Rswitch2)then1140:           if(r2 .gt. Rswitch2)then
1096:            SAT=SA(POS)1141:            SAT=SA(POS)
1097:            SBT=SB(POS)1142:            SBT=SB(POS)
1098:            f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)1143:            f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)
1099:            energy=energy-SAT/3.0*RD1**3-SBT/4.0*RD1**41144:            energy=energy-SAT/3.0*RD1**3-SBT/4.0*RD1**4
1100: !          elseif(r2 .lt. Rswitch2)then1145: !          elseif(r2 .lt. Rswitch2)then
1101:          ! normal repulsive term1146:          ! normal repulsive term
1102: !          else1147: !          else
1103:          ! things should have fallen in one of the previous two...1148:          ! things should have fallen in one of the previous two...
1104: !           write(*,*) 'something went wrong with switching function'1149: !           write(*,*) 'something went wrong with switching function'
1105:           endif1150:           endif
1106:           grad(C13-2) = grad(C13-2) - f_over_r * dx1151:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx
1107:           grad(C13-1) = grad(C13-1) - f_over_r * dy1152:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy
1108:           grad(C13)   = grad(C13)   - f_over_r * dz1153:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz
1109:           grad(C23-2) = grad(C23-2) + f_over_r * dx1154:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx
1110:           grad(C23-1) = grad(C23-1) + f_over_r * dy1155:           grad(C2*3-1) = grad(C2*3-1) + f_over_r * dy
1111:           grad(C23)   = grad(C23)   + f_over_r * dz1156:           grad(C2*3)   = grad(C2*3)   + f_over_r * dz
1112:          endif1157:          endif
1113:       ENDDO1158:       ENDDO
1114: !$OMP END DO1159: !$OMP END DO
1115: !$OMP END PARALLEL1160: !$OMP END PARALLEL
1116:       END1161:       END
1117: 1162: 
1118: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^1163: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
1119: 1164: 
1120: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1165: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1121: ! FUNCTIONS for DH Calculations1166: ! FUNCTIONS for DH Calculations
1133:       implicit NONE1178:       implicit NONE
1134:       double precision kappa,r1179:       double precision kappa,r
1135:       SBMDHVP=-kappa*exp(-kappa*r)/r-exp(-kappa*r)/(r**2)1180:       SBMDHVP=-kappa*exp(-kappa*r)/r-exp(-kappa*r)/(r**2)
1136:       END1181:       END
1137: !****************************************************1182: !****************************************************
1138: 1183: 
1139: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1184: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1140: !* SBMDHELEC computes the forces due to DH electrostatics        *1185: !* SBMDHELEC computes the forces due to DH electrostatics        *
1141: !**********************************************************************1186: !**********************************************************************
1142: 1187: 
1143:       subroutine SBMDHELEC(grad,energy,natoms)1188:       subroutine SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1,NNEXL2,
 1189:      Q NEXCLUSIONS,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
1144: 1190: 
1145:       USE KEY1191:       USE KEY
1146:       USE SBMDATA 
1147:       implicit NONE1192:       implicit NONE
1148:       integer I, N, J, AN, NATOMS1193:       integer I, N, J, AN, NATOMS
1149: 1194: 
1150:       DOUBLE PRECISION grad(3*NATOMS), energy,1195: 
 1196:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
 1197:      Q grad(3*NATOMS), energy,STT,SA,SB,SC,
1151:      Q SBMDHVP, SBMDHV1198:      Q SBMDHVP, SBMDHV
1152:       integer C1,C2,C13,C23,ii,jj,kk,k,l,iii,jjj1199:       integer C1, C2, C12,ii,jj,kk,k,l,iii,jjj
1153:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r 1200:       DOUBLE PRECISION  r1,r2, rm2, rm14, f_over_r, NCswitch,NCcut 
1154:       DOUBLE PRECISION A,B,C, D, COEF2, COEF31201:       DOUBLE PRECISION A,B,C, D, COEF2, COEF3
1155:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,1202:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,
1156:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS1203:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
 1204:         INTEGER NEXCLUSIONS
 1205:         INTEGER NNEXL1(NEXCLUSIONS)
 1206:         INTEGER NNEXL2(NEXCLUSIONS)
1157:         DOUBLE PRECISION dx,dy,dz1207:         DOUBLE PRECISION dx,dy,dz
1158:         integer tempN, alpha1208:         integer tempN, alpha,SBMCHARGEON(NATOMS+1)
1159:         double precision diff,Vfunc,Ffunc1209:         double precision diff,Vfunc,Ffunc,SBMCHARGE(NATOMS)
1160:         double precision DHswitch2,DHcut21210:         double precision PREFACTOR,KAPPA,DHswitch,DHswitch2,DHcut,DHcut2
1161:         ! Ngrid is the number of atoms in that grid point1211:         ! Ngrid is the number of atoms in that grid point
1162:         ! grid is the array of atoms in each grid1212:         ! grid is the array of atoms in each grid
1163:         integer Ngrid,grid,maxgrid,maxpergrid1213:         integer Ngrid,grid,maxgrid,maxpergrid
1164:         ! number of atoms per grid, max1214:         ! number of atoms per grid, max
1165:         parameter (maxpergrid=50)1215:         parameter (maxpergrid=50)
1166:         ! dimensions of grid1216:         ! dimensions of grid
1167:         parameter (maxgrid=15)1217:         parameter (maxgrid=15)
1168:         dimension Ngrid(maxgrid,maxgrid,maxgrid),1218:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
1169:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)1219:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
1170:         integer MaxGridX,MaxGridY,MaxGridZ1220:         integer MaxGridX,MaxGridY,MaxGridZ
1171:         double precision gridsize,RD11221:         double precision gridsize,RD1
1172:         double precision minX,minY,minZ,maxX,maxY,maxZ1222:         double precision minX,minY,minZ,maxX,maxY,maxZ
1173:         integer Xgrid,Ygrid,Zgrid,TID1223:         integer Xgrid,Ygrid,Zgrid,TID
1174:         double precision C2T1224:         double precision C2T
1175: 1225: 
1176:       if(NUMOFSBMCHARGES .eq. 0) RETURN1226:       if(SBMCHARGEON(1) .eq. 1) RETURN
1177: 1227: 
1178:         diff=DHcut-DHswitch1228:         diff=DHcut-DHswitch
1179:         alpha=121229:         alpha=12
1180: 1230: 
1181: 1231: 
1182:         GRIDSIZE=DHCUT*1.011232:         GRIDSIZE=DHCUT*1.01
1183:         DHcut2=DHcut**21233:         DHcut2=DHcut**2
1184:         DHswitch2=DHswitch**21234:         DHswitch2=DHswitch**2
1185: 1235: 
1186:         A=-PREFACTOR*SBMDHV(kappa,DHcut)1236:         A=-PREFACTOR*SBMDHV(kappa,DHcut)
1192: 1242: 
1193: ! grid the system1243: ! grid the system
1194:         minX=100000001244:         minX=10000000
1195:         minY=100000001245:         minY=10000000
1196:         minZ=100000001246:         minZ=10000000
1197:         1247:         
1198:         maxX=-100000001248:         maxX=-10000000
1199:         maxY=-100000001249:         maxY=-10000000
1200:         maxZ=-100000001250:         maxZ=-10000000
1201: 1251: 
1202:         do ii=1,NUMOFSBMCHARGES1252:         do ii=2,SBMCHARGEON(1)
1203:            i=SBMCHARGEON(ii)1253:            i=SBMCHARGEON(ii)
1204:            if(XSBM(i) .lt. minX)then1254: 
1205:                 minX=XSBM(i)1255:            if(X(i) .gt. maxX)then
1206:            elseif(XSBM(i) .gt. maxX)then1256:                 maxX=X(i)
1207:                 maxX=XSBM(i)1257:            elseif(Y(i) .gt. maxY)then
 1258:                 maxY=Y(i)
1208:            endif1259:            endif
1209:            if(YSBM(i) .lt. minY)then1260: 
1210:                 minY=YSBM(i)1261:            if(Z(i) .gt. maxZ)then
1211:            elseif(YSBM(i) .gt. maxY)then1262:                 maxZ=Z(i)
1212:                 maxY=YSBM(i)1263:            elseif(X(i) .lt. minX)then
 1264:                 minX=X(i)
1213:            endif1265:            endif
1214:            if(ZSBM(i) .lt. minZ)then1266: 
1215:                 minZ=ZSBM(i)1267:            if(Y(i) .lt. minY)then
1216:            elseif(ZSBM(i) .gt. maxZ)then1268:                 minY=Y(i)
1217:                 maxZ=ZSBM(i)1269:            elseif(Z(i) .lt. minZ)then
 1270:                 minZ=Z(i)
1218:            endif1271:            endif
1219:         enddo1272:         enddo
1220: 1273: 
1221:         maxgridX=int((maxX-minX)/gridsize)+11274:         maxgridX=int((maxX-minX)/gridsize)+1
1222:         maxgridY=int((maxY-minY)/gridsize)+11275:         maxgridY=int((maxY-minY)/gridsize)+1
1223:         maxgridZ=int((maxZ-minZ)/gridsize)+11276:         maxgridZ=int((maxZ-minZ)/gridsize)+1
1224: 1277: 
1225:         if(maxgridX .ge. maxgrid .or. 1278:         if(maxgridX .ge. maxgrid .or. 
1226:      Q  maxgridY .ge. maxgrid .or.1279:      Q  maxgridY .ge. maxgrid .or.
1227:      Q  maxgridZ .ge. maxgrid )then1280:      Q  maxgridZ .ge. maxgrid )then
1234: 1287: 
1235:         do i=1,maxgrid1288:         do i=1,maxgrid
1236:          do j=1,maxgrid1289:          do j=1,maxgrid
1237:           do k=1,maxgrid1290:           do k=1,maxgrid
1238:                 Ngrid(i,j,k)=01291:                 Ngrid(i,j,k)=0
1239:           enddo1292:           enddo
1240:          enddo1293:          enddo
1241:         enddo1294:         enddo
1242: 1295: 
1243: 1296: 
1244:         do ii=1,NUMOFSBMCHARGES1297:         do ii=2,SBMCHARGEON(1)
1245:                 i=SBMCHARGEON(ii)1298:                 i=SBMCHARGEON(ii)
1246: 1299: 
1247:                 Xgrid=int((XSBM(i)-minX)/gridsize)+11300:                 Xgrid=int((X(i)-minX)/gridsize)+1
1248:                 Ygrid=int((YSBM(i)-minY)/gridsize)+11301:                 Ygrid=int((Y(i)-minY)/gridsize)+1
1249:                 Zgrid=int((ZSBM(i)-minZ)/gridsize)+11302:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
1250:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+11303:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
1251:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then1304:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
1252:                         write(*,*) 'ERROR: too many atoms in a grid in DH'1305:                         write(*,*) 'ERROR: too many atoms in a grid in DH'
1253:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid1306:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
1254:                         STOP1307:                         STOP
1255:                 endif1308:                 endif
1256:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i1309:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1257:         enddo1310:         enddo
1258: 1311: 
1259:         tempN=01312: 
 1313:            tempN=0
 1314: 
1260: 1315: 
1261: !$OMP PARALLEL1316: !$OMP PARALLEL
1262: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,C13,C23,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  1317: !$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)  
1263: !$OMP& REDUCTION(+:energy,grad)1318: !$OMP& REDUCTION(+:energy,grad)
1264:         TID=01319:         TID=0
1265:         NTHREADS=1;1320:         NTHREADS=1;
1266: !$      TID = OMP_GET_THREAD_NUM()1321: !$      TID = OMP_GET_THREAD_NUM()
1267: !$      NTHREADS = OMP_GET_NUM_THREADS()1322: !$      NTHREADS = OMP_GET_NUM_THREADS()
1268:          do iii=1+TID,NUMOFSBMCHARGES,NTHREADS1323:          do iii=2+TID,SBMCHARGEON(1),NTHREADS
1269:           C1=SBMCHARGEON(iii)1324:           C1=SBMCHARGEON(iii)
1270:           C13=C1*31325:           Xgrid=int((X(C1)-minX)/gridsize)+1
1271:           Xgrid=int((XSBM(C1)-minX)/gridsize)+11326:           Ygrid=int((Y(C1)-minY)/gridsize)+1
1272:           Ygrid=int((YSBM(C1)-minY)/gridsize)+11327:           Zgrid=int((Z(C1)-minZ)/gridsize)+1
1273:           Zgrid=int((ZSBM(C1)-minZ)/gridsize)+1 
1274:         1328:         
1275:           do ii=XGRID-1,XGRID+11329:           do ii=XGRID-1,XGRID+1
1276:            do jj=YGRID-1,YGRID+11330:            do jj=YGRID-1,YGRID+1
1277:             do kk=ZGRID-1,ZGRID+11331:             do kk=ZGRID-1,ZGRID+1
1278:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1332:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1279:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1333:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1280:               do jjj=1,Ngrid(ii,jj,kk)1334:               do jjj=1,Ngrid(ii,jj,kk)
1281:                C2=grid(ii,jj,kk,jjj)1335:                C2=grid(ii,jj,kk,jjj)
1282:                C23=C2*31336:            !C12=C2-C1
 1337:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNEXL(C1,C12))) )then
1283:                if(C2-C1 .gt. 0)then1338:                if(C2-C1 .gt. 0)then
1284:                 dx = XSBM(C1) - XSBM(C2)1339:                 dx = X(C1) - X(C2)
1285:                 dy = YSBM(C1) - YSBM(C2)1340:                 dy = Y(C1) - Y(C2)
1286:                 dz = ZSBM(C1) - ZSBM(C2)1341:                 dz = Z(C1) - Z(C2)
1287:                 r2 = dx**2 + dy**2 + dz**21342:                 r2 = dx**2 + dy**2 + dz**2
1288:                 if(r2 .le. DHcut2)then1343:                 if(r2 .le. DHcut2)then
1289:                  C2T=SBMCHARGE(C1)*SBMCHARGE(C2)1344:                  C2T=SBMCHARGE(C1)*SBMCHARGE(C2)
1290:         ! add force evaluation now1345:         ! add force evaluation now
1291:                  rm2 = 1/r21346:                  rm2 = 1/r2
1292:                  r1=sqrt(r2)1347:                  r1=sqrt(r2)
1293:                  energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)1348:                  energy=energy+PREFACTOR*C2T*SBMDHV(kappa,r1)
1294:                  f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)1349:                  f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)
1295:                  if(r2 .gt. DHswitch2)then1350:                  if(r2 .gt. DHswitch2)then
1296:                   RD1=r1-DHswitch1351:                   RD1=r1-DHswitch
1297:                   f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))1352:                   f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))
1298:                   energy=energy+COEF2*RD1**2+COEF3*RD1**31353:                   energy=energy+COEF2*RD1**2+COEF3*RD1**3
1299:                  elseif(r2 .lt. DHswitch2)then1354:                  elseif(r2 .lt. DHswitch2)then
1300:                 ! normal DH term1355:                 ! normal DH term
1301:                  else1356:                  else
1302:                 ! things should have fallen in one of the previous two...1357:                 ! things should have fallen in one of the previous two...
1303:                   write(*,*) 'something went wrong with DH switching function'1358:                   write(*,*) 'something went wrong with DH switching function'
1304:                   STOP1359:                   STOP
1305:                  endif1360:                  endif
1306:                  f_over_r=f_over_r*sqrt(rm2)1361:                  f_over_r=f_over_r*sqrt(rm2)
1307:                  grad(C13-2) = grad(C13-2) + f_over_r * dx1362:                  grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx
1308:                  grad(C13-1) = grad(C13-1) + f_over_r * dy1363:                  grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy
1309:                  grad(C13)   = grad(C13)   + f_over_r * dz1364:                  grad(C1*3)   = grad(C1*3)   + f_over_r * dz
1310:                  grad(C23-2) = grad(C23-2) - f_over_r * dx1365:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
1311:                  grad(C23-1) = grad(C23-1) - f_over_r * dy1366:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
1312:                  grad(C23)   = grad(C23)   - f_over_r * dz1367:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz
1313: 1368: 
1314:                   endif1369:                   endif
1315:                 endif1370:                 endif
1316:                enddo1371:                enddo
1317:                   endif1372:                   endif
1318:              enddo1373:              enddo
1319:             enddo1374:             enddo
1320:            enddo1375:            enddo
1321:           enddo1376:           enddo
1322: 1377: 
1323: !!! Add routine to subtract exclusion interactions1378: !!! Add routine to subtract exclusion interactions
1324: !$OMP DO1379: !$OMP DO
1325:        DO I=1,NEXCLUSIONS1380:        DO I=1,NEXCLUSIONS
1326:          C1=NNEXL1(I) 1381:          C1=NNEXL1(I) 
1327:          C2=NNEXL2(I)1382:          C2=NNEXL2(I)
1328:          C13=C1*31383:          dx = X(C1) - X(C2)
1329:          C23=C2*31384:          dy = Y(C1) - Y(C2)
1330:          dx = XSBM(C1) - XSBM(C2)1385:          dz = Z(C1) - Z(C2)
1331:          dy = YSBM(C1) - YSBM(C2) 
1332:          dz = ZSBM(C1) - ZSBM(C2) 
1333:          r2 = dx**2 + dy**2 + dz**21386:          r2 = dx**2 + dy**2 + dz**2
1334:          if(r2 .le. DHcut2)then1387:          if(r2 .le. DHcut2)then
1335:           C2T=SBMCHARGE(C1)*SBMCHARGE(C2)1388:           C2T=SBMCHARGE(C1)*SBMCHARGE(C2)
1336:           rm2 = 1/r21389:           rm2 = 1/r2
1337:           r1=sqrt(r2)1390:           r1=sqrt(r2)
1338:           energy=energy-PREFACTOR*C2T*SBMDHV(kappa,r1)1391:           energy=energy-PREFACTOR*C2T*SBMDHV(kappa,r1)
1339:           f_over_r=-PREFACTOR*C2T*SBMDHVP(kappa,r1)1392:           f_over_r=-PREFACTOR*C2T*SBMDHVP(kappa,r1)
1340:           if(r2 .gt. DHswitch2)then1393:           if(r2 .gt. DHswitch2)then
1341:            RD1=r1-DHswitch1394:            RD1=r1-DHswitch
1342:            f_over_r=f_over_r-C2T*(2*COEF2*RD1+3*COEF3*RD1**2)1395:            f_over_r=f_over_r-C2T*(2*COEF2*RD1+3*COEF3*RD1**2)
1343:            energy=energy-COEF2*RD1**2+COEF3*RD1**31396:            energy=energy-COEF2*RD1**2+COEF3*RD1**3
1344:           elseif(r2 .lt. DHswitch2)then1397:           elseif(r2 .lt. DHswitch2)then
1345:          ! normal DH term1398:          ! normal DH term
1346:           else1399:           else
1347:          ! things should have fallen in one of the previous two...1400:          ! things should have fallen in one of the previous two...
1348:            write(*,*) 'something went wrong with DH switching function'1401:            write(*,*) 'something went wrong with DH switching function'
1349:           endif1402:           endif
1350:           f_over_r=f_over_r*sqrt(rm2)1403:           f_over_r=f_over_r*sqrt(rm2)
1351:           grad(C13-2) = grad(C13-2) + f_over_r * dx1404:           grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx
1352:           grad(C13-1) = grad(C13-1) + f_over_r * dy1405:           grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy
1353:           grad(C13)   = grad(C13)   + f_over_r * dz1406:           grad(C1*3)   = grad(C1*3)   + f_over_r * dz
1354:           grad(C23-2) = grad(C23-2) - f_over_r * dx1407:           grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
1355:           grad(C23-1) = grad(C23-1) - f_over_r * dy1408:           grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
1356:           grad(C23)   = grad(C23)   - f_over_r * dz1409:           grad(C2*3)   = grad(C2*3)   - f_over_r * dz
1357:          ENDIF1410:          ENDIF
1358:        ENDDO1411:        ENDDO
1359: !$OMP END DO1412: !$OMP END DO
1360: !$OMP END PARALLEL1413: !$OMP END PARALLEL
1361: 1414: 
 1415: 
 1416: 
1362:       END1417:       END
1363: 1418: 
1364: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^1419: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMDHELEC^^^^^^^^^^^^^^^^^^^^^^^^^
1365: 1420: 
1366: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1421: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1367: !* SBMPR computes the forces due to position restraints               *1422: !* SBMPR computes the forces due to position restraints               *
1368: !**********************************************************************1423: !**********************************************************************
1369: 1424: 
1370:       subroutine SBMPR(grad, energy, natoms)1425:       subroutine SBMPR(x,y,z,grad, energy, natoms, 
 1426:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
1371: 1427: 
1372:       USE KEY1428:       USE KEY
1373:       USE SBMDATA 
1374:       implicit NONE1429:       implicit NONE
1375:       integer I, J, J3,NATOMS,I1,I21430:       integer I, J, NATOMS,I1,I2
1376:       DOUBLE PRECISION grad(3*NATOMS), energy,ET1431:       INTEGER SBMPRN, SBMPRI(NATOMS)
 1432:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)
 1433:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
 1434:      Q grad(3*NATOMS), energy,ET
1377:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K61435:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K6
1378: 1436: 
1379: 1437: 
1380:       if(SBMPRN .eq. 0) RETURN1438:       if(SBMPRN .eq. 0) RETURN
1381: 1439: 
1382: !$OMP PARALLEL PRIVATE(I,J,J3,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:energy,grad)1440: !$OMP PARALLEL PRIVATE(I,J,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:energy,grad)
1383: !$OMP DO1441: !$OMP DO
1384:               do I=1,SBMPRN1442:               do I=1,SBMPRN
1385:            J=SBMPRI(I)1443:            J=SBMPRI(I)
1386:            I1=(I-1)*61444:            I1=(I-1)*6
1387:            I2=(I-1)*31445:            I2=(I-1)*3
1388:            J3=J*31446:            DX=X(J)-SBMPRX(I2+1)
1389:            DX=XSBM(J)-SBMPRX(I2+1)1447:            DY=Y(J)-SBMPRX(I2+2)
1390:            DY=YSBM(J)-SBMPRX(I2+2)1448:            DZ=Z(J)-SBMPRX(I2+3)
1391:            DZ=ZSBM(J)-SBMPRX(I2+3) 
1392:            K1=SBMPRK(I1+1)1449:            K1=SBMPRK(I1+1)
1393:            K2=SBMPRK(I1+2)1450:            K2=SBMPRK(I1+2)
1394:            K3=SBMPRK(I1+3)1451:            K3=SBMPRK(I1+3)
1395:            K4=SBMPRK(I1+4)1452:            K4=SBMPRK(I1+4)
1396:            K5=SBMPRK(I1+5)1453:            K5=SBMPRK(I1+5)
1397:            K6=SBMPRK(I1+6)1454:            K6=SBMPRK(I1+6)
1398:            1455:            
1399:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +1456:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +
1400:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ1457:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ
1401: 1458: 
1402:            grad(J3-2) = grad(J3-2) + K1*DX + K4*DY + K5*DZ 1459:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ 
1403:            grad(J3-1) = grad(J3-1) + K2*DY + K4*DX + K6*DZ1460:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ
1404:            grad(J3)   = grad(J3)   + K3*DZ + K5*DX + K6*DY1461:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY
1405:         enddo1462:         enddo
1406: !$OMP END DO1463: !$OMP END DO
1407: !$OMP END PARALLEL1464: !$OMP END PARALLEL
1408: 1465: 
1409:       END1466:       END
1410: 1467: 
1411: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^1468: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMPR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1412: 1469: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0