hdiff output

r30688/SBM.f 2016-07-06 15:35:57.071515506 +0100 r30687/SBM.f 2016-07-06 15:35:57.435520430 +0100
  3:       implicit NONE  3:       implicit NONE
  4:       INTEGER NATOMS,i,j  4:       INTEGER NATOMS,i,j
  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)
  6:       DOUBLE PRECISION ENERGY,STT  6:       DOUBLE PRECISION ENERGY,STT
  7:   7: 
  8:       LOGICAL :: CALLED=.FALSE.  8:       LOGICAL :: CALLED=.FALSE.
  9:       LOGICAL GTEST, STEST  9:       LOGICAL GTEST, STEST
 10:         integer NSBMMAX 10:         integer NSBMMAX
 11:         parameter(NSBMMAX=30000) 11:         parameter(NSBMMAX=30000)
 12:       INTEGER MAXSEP,MAXSEPSYS 12:       INTEGER MAXSEP,MAXSEPSYS
 13:       PARAMETER (MAXSEP=40) 13:       PARAMETER (MAXSEP=20)
 14:       INTEGER MAXATOMTYPES 
 15:       INTEGER ATOMTYPES 
 16:       PARAMETER (MAXATOMTYPES=10) 
 17:       DIMENSION STT(MAXATOMTYPES*MAXATOMTYPES),ATOMTYPES(NSBMMAX) 
 18:       DOUBLE PRECISION SA, SB, SC 
 19:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES) 
 20:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES) 
 21:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES) 
 22:       LOGICAL NNCINC 14:       LOGICAL NNCINC
 23:       DIMENSION NNCINC(NSBMMAX*MAXSEP) 15:       DIMENSION NNCINC(NSBMMAX*MAXSEP)
 24:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 16:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX),
 25:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF(NSBMMAX*5*6), 17:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF(NSBMMAX*5*6),
 26:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 18:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
 27:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2), 19:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2),
 28:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3), 20:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3),
 29:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5), 21:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5),
 30:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX) 22:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX)
 31:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NEXCLUDE,NEXCLUDEELEC 23:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NEXCLUDE,NEXCLUDEELEC
 32:       PARAMETER (MAXEXCLUSIONS=20) 24:       PARAMETER (MAXEXCLUSIONS=20)
 33:       PARAMETER (MAXEXCLUSIONSELEC=5) 25:       PARAMETER (MAXEXCLUSIONSELEC=5)
 34:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS), 26:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS),
 35:      Q NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC) 27:      Q NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC)
 36:  28: 
 37:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3) 29:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3)
 38:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX) 30:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX)
 39:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 31:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 40:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 32:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
 41:      Q SBMPRK,SBMPRX,CONCOEF,SA,SB,SC 33:      Q SBMPRK,SBMPRX,CONCOEF
 42:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC, 34:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC,
 43:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1, 35:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,
 44:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,MAXSEPSYS, 36:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,MAXSEPSYS
 45:      Q ATOMTYPES 
 46:         common /logical/ NNCINC 37:         common /logical/ NNCINC
 47:  38: 
 48:         if(NATOMS.gt. NSBMMAX)then 39:         if(NATOMS.gt. NSBMMAX)then
 49:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 40:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 50:         STOP 41:         STOP
 51:         endif 42:         endif
 52:  43: 
 53:        if(.NOT.CALLED)then 44:        if(.NOT.CALLED)then
 54:         NEXCLUDE=MAXEXCLUSIONS 45:         NEXCLUDE=MAXEXCLUSIONS
 55:         NEXCLUDEELEC=MAXEXCLUSIONS 46:         NEXCLUDEELEC=MAXEXCLUSIONS
 56:         MAXSEPSYS=0 47:         MAXSEPSYS=0
 57:         write(*,*) 48:         write(*,*)
 58:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 49:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 59:         write(*,*)  '   For more information on SMOG models, see:' 50:         write(*,*)  '   For more information on SMOG models, see:'
 60:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 51:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 61:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 52:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 62:         write(*,*) 53:         write(*,*)
 63:          54:         
 64:  55: 
 65:        call SBMinit(NATOMS,ATOMTYPES,MAXATOMTYPES,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP, 56:        call SBMinit(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,
 66:      Q PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE, 57:      Q PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE,
 67:      Q NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC, 58:      Q NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,
 68:      Q SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,SA,SB,SC,CONTACTTYPE, 59:      Q SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,CONTACTTYPE,
 69:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 60:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 70:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or. 61:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or.
 71:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then 62:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then
 72:         write(*,*) 'ERROR: arrays not large large enough for the system' 63:         write(*,*) 'ERROR: arrays not large large enough for the system'
 73:         write(*,*) '   entries=',NBA, NTA,NPA,NC  64:         write(*,*) '   entries=',NBA, NTA,NPA,NC 
 74:         stop 65:         stop
 75:         endif  66:         endif 
 76:         CALLED=.TRUE. 67:         CALLED=.TRUE.
 77:         endIF 68:         endIF
 78:  69: 
 79:  70: 
 80: ! call the energy routine 71: ! call the energy routine
 81:       call calc_energy_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,energy, 72:       call calc_energy_SBM(qo,natoms,GRAD,energy,Ib1,Ib2,Rb,bK,IT,JT,KT,
 82:      Q Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC, 73:      Q ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF, 
 83:      Q JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC, 74:      Q NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,
 84:      Q NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE, 75:      Q MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,STT,NCcut,
 85:      Q SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,KAPPA, 76:      Q CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,
 86:      Q DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 77:      Q SBMPRX)
 87:       IF (STEST) THEN 78:       IF (STEST) THEN
 88:          PRINT '(A)','ERROR - second derivatives not available' 79:          PRINT '(A)','ERROR - second derivatives not available'
 89:          STOP 80:          STOP
 90:       ENDIF 81:       ENDIF
 91:       return 82:       return
 92:       end 83:       end
 93:  84: 
 94:  85: 
 95: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 86: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 96: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 87: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 97: !* startt then the velocities are assigned, otherwise, they are read   * 88: !* startt then the velocities are assigned, otherwise, they are read   *
 98: !* by selecting 2, or generated by selecting 3                         * 89: !* by selecting 2, or generated by selecting 3                         *
 99: !*********************************************************************** 90: !***********************************************************************
100:  91: 
101:       subroutine SBMINIT(NATOMS,ATOMTYPES,MAXATOMTYPES,Ib1,Ib2,Rb,bK,IT, 92:       subroutine SBMINIT(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,
102:      Q JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF,  93:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF, 
103:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNEXL1ELEC,NNEXL2ELEC,MAXEXCLUSIONSELEC, 94:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNEXL1ELEC,NNEXL2ELEC,MAXEXCLUSIONSELEC,
104:      Q NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON, 95:      Q NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,
105:      Q NCswitch,NCcut,STT,SA,SB,SC,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch, 96:      Q NCswitch,NCcut,STT,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,
106:      Q DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 97:      Q DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
107:       USE KEY 98:       USE KEY
108:       USE COMMONS, only: ATMASS 99:       USE COMMONS, only: ATMASS
109:       implicit NONE100:       implicit NONE
110: 101: 
111:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,102:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,
112:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 103:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 
113:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,104:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,
114:      Q MAXEXCLUSIONS,MAXEXCLUSIONSNEW,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,105:      Q MAXEXCLUSIONS,MAXEXCLUSIONSNEW,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,MAXSEPSYS
115:      Q MAXSEPSYS,ATOMTYPES(NATOMS),ATYPE,POS,MAXATOMTYPES,NATOMTYPES106: 
116:       DOUBLE PRECISION :: Rb(2*NATOMS), bK(2*NATOMS), ANTC(2*NATOMS), 107:       DOUBLE PRECISION :: Rb(2*NATOMS), bK(2*NATOMS), ANTC(2*NATOMS), 
117:      Q Tk(2*NATOMS), PK(3*NATOMS), PHISBM(3*NATOMS),108:      Q Tk(2*NATOMS), PK(3*NATOMS), PHISBM(3*NATOMS),
118:      Q Sigma, EpsC, NCswitch,NCcut,TMPSTT(MAXATOMTYPES),109:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS)
119:      Q STT(MAXATOMTYPES**2),SBMCHARGE(NATOMS)110: !      DOUBLE PRECISION CONCOEF
120:       DOUBLE PRECISION RDIFF,ALPHA111:       DOUBLE PRECISION CONCOEF (NC*6)
121:       DOUBLE PRECISION CONCOEF(NC*6) 
122:       INTEGER  Ib1(2*NATOMS), Ib2(2*NATOMS), IT(2*NATOMS), JT(2*NATOMS), 112:       INTEGER  Ib1(2*NATOMS), Ib2(2*NATOMS), IT(2*NATOMS), JT(2*NATOMS), 
123:      Q KT(2*NATOMS),IP(3*NATOMS),JP(3*NATOMS),KP(3*NATOMS),113:      Q KT(2*NATOMS),IP(3*NATOMS),JP(3*NATOMS),KP(3*NATOMS),
124:      Q LP(3*NATOMS),IC(NATOMS*5),JC(NATOMS*5), 114:      Q LP(3*NATOMS),IC(NATOMS*5),JC(NATOMS*5), 
125:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE115:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE
126:       INTEGER PHITYPE(3*NATOMS),M116:       INTEGER PHITYPE(3*NATOMS),M
127:       integer AA,BB,ANTEMP117:       integer AA,BB,ANTEMP
128:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)118:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)
129:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)119:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)
130:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)120:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)
131:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)121:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)
132:       DOUBLE PRECISION dx,dy,dz122:       DOUBLE PRECISION dx,dy,dz
133:       double precision PI,ARG1,ARG2,ARG3,ARG4123:       double precision PI,ARG1,ARG2,ARG3,ARG4
134:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut124:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
135:       integer CONTACTTYPE(NATOMS*5),TMPARG125:       integer CONTACTTYPE(NATOMS*5),TMPARG
136:       DOUBLE PRECISION SA, SB, SC126:  
137:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES) 
138:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES) 
139:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES) 
140:  
141:       character TMPCHAR127:       character TMPCHAR
142:       integer TMPINT,NUMCHARGES,nexc,I1,I2128:       integer TMPINT,NUMCHARGES,nexc,I1,I2
143:       double precision TMPREAL,concentration129:       double precision TMPREAL,concentration
144:       LOGICAL NNCINC(NATOMS*MAXSEP),TARR(NATOMS)130:       LOGICAL NNCINC(NATOMS*MAXSEP),TARR(NATOMS)
145: 131: 
146:       INTEGER SBMPRN, SBMPRI(NATOMS)132:       INTEGER SBMPRN, SBMPRI(NATOMS)
147:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)133:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)
148:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS134:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS
149:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)135:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)
150: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 136: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
170: 156: 
171: ! These lines read in the parameters.157: ! These lines read in the parameters.
172:         open(30, file='SBM.INP', status='old', access='sequential')158:         open(30, file='SBM.INP', status='old', access='sequential')
173:         write(*,*) 159:         write(*,*) 
174:         write(*,*) 'Reading the forcefield from SBM.INP'160:         write(*,*) 'Reading the forcefield from SBM.INP'
175:         write(*,*) 161:         write(*,*) 
176:         read(30,*)162:         read(30,*)
177:         read(30,*)163:         read(30,*)
178:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut164:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut
179:         read(30,*)165:         read(30,*)
180:         read(30,*) NCswitch, NCcut166:         read(30,*) RSig, Reps, NCswitch, NCcut
181:         write(*,*) 'NCswitch,NCcut'167:         write(*,*) 'RSig, Reps,NCswitch,NCcut'
182:         write(*,'(2F10.5)') NCswitch, NCcut168:         write(*,'(4F10.5)') RSig, Reps, NCswitch, NCcut
183:         write(*,*) 'Reading electrostatic parameters'169:         write(*,*) 'Reading electrostatic parameters'
184:         write(*,'(5F10.5)') PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut170:         write(*,'(5F10.5)') PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut
185:         read(30,*) NATOMTYPES 
186:         do I=1,NATOMTYPES 
187:           read(30,*) TMPARG,rsig,reps 
188:           if(TMPARG .ne. I)then 
189:             write(*,*) 'ERROR: atomtypes must be sequential integers' 
190:             STOP   
191:           endif 
192:           TMPSTT(TMPARG)=reps*rsig**12 
193:         enddo 
194:         ! using combination rule 1 
195:  
196:         Rdiff=NCcut-NCswitch 
197:         alpha=12 
198:  
199:         do I=1,NATOMTYPES 
200:           do J=1,NATOMTYPES 
201:             POS=(I-1)*MAXATOMTYPES+J 
202:             STT(POS)=sqrt(TMPSTT(I)*TMPSTT(J)) 
203:             SB(POS)=-1.0/Rdiff**3*( 2*alpha*STT(POS)/NCcut**(alpha+1)  +  
204:      Q         (alpha)*(alpha+1)*STT(POS)*Rdiff/NCcut**(alpha+2)) 
205:             SA(POS)=-(alpha*(alpha+1)*STT(POS)/NCcut**(alpha+2)+ 
206:      Q         3*SB(POS)*Rdiff**2)/(2*Rdiff) 
207:             SC(POS)=-(STT(POS)/NCcut**alpha + 
208:      Q         SA(POS)/3.0*Rdiff**3+SB(POS)/4.0*Rdiff**4) 
209:           enddo  
210:         enddo 
211:  
212:  
213:  
214:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is in units171:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is in units
215:         ! A^-1172:         ! A^-1
216:         KAPPA=0.32*sqrt(CONCENTRATION)173:         KAPPA=0.32*sqrt(CONCENTRATION)
217:         PREFACTOR=PREFACTOR/DC174:         PREFACTOR=PREFACTOR/DC
218:         read(30,*) ANtemp175:         read(30,*) ANtemp
219:         write(*,*) ANtemp, 'atoms'176:         write(*,*) ANtemp, 'atoms'
220:         NUMCHARGES=1177:         NUMCHARGES=1
221:         do i=1, ANtemp178:         do i=1, ANtemp
222:           read(30,*) TMPINT,ATYPE,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)179:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)
223:             if(ATYPE .gt. NATOMTYPES)THEN 
224:               write(*,*) 'ERROR: Unknown atomtype',ATYPE 
225:             endif 
226:             ATOMTYPES(I)=ATYPE 
227:           if(SBMCHARGE(i) .ne. 0)then180:           if(SBMCHARGE(i) .ne. 0)then
228:              NUMCHARGES=NUMCHARGES+1181:              NUMCHARGES=NUMCHARGES+1
229:              SBMCHARGEON(NUMCHARGES)=i182:              SBMCHARGEON(NUMCHARGES)=i
230:           endif 183:           endif 
231:         end do184:         end do
232:         SBMCHARGEON(1)=NUMCHARGES185:         SBMCHARGEON(1)=NUMCHARGES
233:         186:         
234:         read(30,*) 187:         read(30,*) 
235:         read(30,*) NC188:         read(30,*) NC
236:         write(*,*) NC, 'contacts'        189:         write(*,*) NC, 'contacts'        
298:           read(30,*) nTA251:           read(30,*) nTA
299:           write(*,*) nTA, 'angles'252:           write(*,*) nTA, 'angles'
300:         do i=1, nTA253:         do i=1, nTA
301:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)254:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)
302:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          255:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
303:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          256:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
304:           CALL INCLUDEEXCLUSIONS(NATOMS,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          257:           CALL INCLUDEEXCLUSIONS(NATOMS,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
305:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.258:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.
306:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.259:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.
307:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN260:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN
308:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS GREATER THAN MAXSEP'261:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS LESS THAN MAXSEP'
309:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'262:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
310:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'263:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
311:           ENDIF264:           ENDIF
312:         ENDDO265:         ENDDO
313: 266: 
314:           read(30,*) 267:           read(30,*) 
315:           read(30,*) nPA268:           read(30,*) nPA
316:           write(*,*) nPA, 'dihedrals'269:           write(*,*) nPA, 'dihedrals'
317: 270: 
318: ! this reads in the dihedral angles and calculates the cosines and sines271: ! this reads in the dihedral angles and calculates the cosines and sines
325:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          278:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
326:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          279:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
327:           CALL INCLUDEEXCLUSIONS(NATOMS,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          280:           CALL INCLUDEEXCLUSIONS(NATOMS,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          
328: 281: 
329:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.282:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.
330:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.283:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.
331:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.284:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.
332:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.285:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.
333:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.286:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.
334:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN287:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN
335:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS GREATER THAN MAXSEP'288:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS LESS THAN MAXSEP'
336:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'289:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
337:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'290:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
338:           ENDIF291:           ENDIF
339:         ENDDO 292:         ENDDO 
340: 293: 
341:         read(30,*)294:         read(30,*)
342:         read(30,*) nexc295:         read(30,*) nexc
343:         write(*,*) nexc, 'exculusions'296:         write(*,*) nexc, 'exculusions'
344: 297: 
345:       do i=1, nexc298:       do i=1, nexc
365:             IF(SBMCHARGE(I) .NE. 0 .AND. SBMCHARGE(M) .NE. 0)THEN318:             IF(SBMCHARGE(I) .NE. 0 .AND. SBMCHARGE(M) .NE. 0)THEN
366:               MAXEXCLUSIONSELEC=MAXEXCLUSIONSELEC+1319:               MAXEXCLUSIONSELEC=MAXEXCLUSIONSELEC+1
367:               IF(MAXEXCLUSIONSELEC .GT. NTE)THEN320:               IF(MAXEXCLUSIONSELEC .GT. NTE)THEN
368:                 write(*,*) 'ERROR: TOO MANY EXCLUSIONS (ELEC), IN TOTAL.'321:                 write(*,*) 'ERROR: TOO MANY EXCLUSIONS (ELEC), IN TOTAL.'
369:               ENDIF322:               ENDIF
370:               NNEXL1ELEC(MAXEXCLUSIONSELEC)=I323:               NNEXL1ELEC(MAXEXCLUSIONSELEC)=I
371:               NNEXL2ELEC(MAXEXCLUSIONSELEC)=M324:               NNEXL2ELEC(MAXEXCLUSIONSELEC)=M
372:             ENDIF325:             ENDIF
373:           ENDDO326:           ENDDO
374:         ENDDO327:         ENDDO
 328:         STT=reps*rsig**12
375: ! At this point, we will change the meaning of maxexclusions, which is329: ! At this point, we will change the meaning of maxexclusions, which is
376: ! actually NEXCLUDE330: ! actually NEXCLUDE
377:         MAXEXCLUSIONS=MAXEXCLUSIONSNEW331:         MAXEXCLUSIONS=MAXEXCLUSIONSNEW
378:         do I=1,NATOMS332:         do I=1,NATOMS
379:           TARR(I)=.FALSE.333:           TARR(I)=.FALSE.
380:         enddo334:         enddo
381: 335: 
382: ! read in position restraints336: ! read in position restraints
383:        read(30,*) 337:        read(30,*) 
384:        read(30,*) SBMPRN338:        read(30,*) SBMPRN
448:           EXCLIST((M-1)*MAXEXCLUDE+NEXCLUDE(M))=N402:           EXCLIST((M-1)*MAXEXCLUDE+NEXCLUDE(M))=N
449:         ENDIF403:         ENDIF
450:       ENDIF404:       ENDIF
451:       END405:       END
452: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^406: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
453: 407: 
454: 408: 
455: C409: C
456: C Calculate the Forces and energies410: C Calculate the Forces and energies
457: C411: C
458:       subroutine CALC_ENERGY_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,412:       subroutine CALC_ENERGY_SBM(qo,natoms,GRAD,energy,Ib1,Ib2,
459:      Q energy,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,413:      Q Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,
460:      Q PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,414:      Q CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,NNEXL2ELEC,
461:      Q NNEXL2ELEC,NEXCLUSIONSELEC,NCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,415:      Q NEXCLUSIONSELEC,NCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,
462:      Q SBMCHARGE,SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,416:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
463:      Q KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)417:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
464: 418: 
465:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,419:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,
466:      Q NEXCLUSIONSELEC,ATOMTYPES(NATOMS),MAXATOMTYPES420:      Q NEXCLUSIONSELEC
467: 421: 
468:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY422:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
469:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)423:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
470:       LOGICAL NNCINC424:       LOGICAL NNCINC
471:       DIMENSION NNCINC(NATOMS*MAXSEP)425:       DIMENSION NNCINC(NATOMS*MAXSEP)
472:       DOUBLE PRECISION SA, SB, SC 
473:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES) 
474:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES) 
475:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES) 
476: 426: 
477:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 427:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
478:      Q PHISBM(NPA),NCswitch,NCcut,STT(MAXATOMTYPES**2), SBMCHARGE(NATOMS)428:      Q PHISBM(NPA),NCswitch,NCcut,STT, 
 429:      Q SBMCHARGE(NATOMS)
479:         DOUBLE PRECISION  CONCOEF(NC*6)430:         DOUBLE PRECISION  CONCOEF(NC*6)
480:         INTEGER Ib1(NBA),Ib2(NBA),IT(NTA),JT(NTA),KT(NTA),IP(NPA), 431:         INTEGER Ib1(NBA),Ib2(NBA),IT(NTA),JT(NTA),KT(NTA),IP(NPA), 
481:      Q JP(NPA),KP(NPA),LP(NPA),IC(NC),JC(NC),SBMCHARGEON(NATOMS),432:      Q JP(NPA),KP(NPA),LP(NPA),IC(NC),JC(NC),SBMCHARGEON(NATOMS),
482:      Q PHITYPE(NPA),CONTACTTYPE(NC)433:      Q PHITYPE(NPA),CONTACTTYPE(NC)
483:         INTEGER NNEXL1(NEXCLUSIONS)434:         INTEGER NNEXL1(NEXCLUSIONS)
484:         INTEGER NNEXL2(NEXCLUSIONS)435:         INTEGER NNEXL2(NEXCLUSIONS)
485:         INTEGER NNEXL1ELEC(NEXCLUSIONSELEC)436:         INTEGER NNEXL1ELEC(NEXCLUSIONSELEC)
486:         INTEGER NNEXL2ELEC(NEXCLUSIONSELEC)437:         INTEGER NNEXL2ELEC(NEXCLUSIONSELEC)
487:       DOUBLE PRECISION dx,dy,dz438:       DOUBLE PRECISION dx,dy,dz
488:       INTEGER SBMPRN, SBMPRI(NATOMS)439:       INTEGER SBMPRN, SBMPRI(NATOMS)
489:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)440:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)
490:       do i = 1, natoms441:       do i = 1, natoms
491:         j = (i-1)*3442:         j = (i-1)*3
492:         x(i) = qo(j+1)443:         x(i) = qo(j+1)
493:         y(i) = qo(j+2)444:         y(i) = qo(j+2)
494:         z(i) = qo(j+3)445:         z(i) = qo(j+3)
495:         grad(j+1) = 0.0446:         grad(j+1) = 0.0
496:         grad(j+2) = 0.0447:         grad(j+2) = 0.0
497:         grad(j+3) = 0.0448:         grad(j+3) = 0.0
498:       enddo449:       enddo
 450: 
499:       energy = 0.0451:       energy = 0.0
500:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)452:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
501:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)453:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
502:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,454:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
503:      Q PHITYPE,PHISBM,NPA)455:      Q PHITYPE,PHISBM,NPA)
504:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC,CONCOEF,456:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC,CONCOEF,
505:      Q NC,CONTACTTYPE)457:      Q NC,CONTACTTYPE)
506:       call SBMNonContacts(x,y,z,grad,energy,natoms,ATOMTYPES,458:       call SBMNonContacts(x,y,z,grad,energy,natoms,NNEXL1,NNEXL2,
507:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,459:      Q NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)
508:      Q NCswitch,NCcut,STT,SA,SB,SC) 
509:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,460:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,
510:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)461:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
511:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)462:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
512: 463: 
513:       end464:       end
514: 465: 
515: 466: 
516: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<467: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
517: !* SBMBonds  computes the hookean force and energy between chosen atoms *468: !* SBMBonds  computes the hookean force and energy between chosen atoms *
518: !***********************************************************************469: !***********************************************************************
922: !$OMP END DO873: !$OMP END DO
923: !$OMP END PARALLEL874: !$OMP END PARALLEL
924: 875: 
925: 876: 
926:       end877:       end
927: 878: 
928: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^879: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
929: 880: 
930: 881: 
931: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<882: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
932: !* SBMNonContacts computes the forces due to non native contacts      *883: !* SBMNonContacts computes the forces due to non native contacts       *
933: !**********************************************************************884: !**********************************************************************
934: 885: 
935:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,ATOMTYPES,886:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,NNEXL1,NNEXL2,
936:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,887:      Q  NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)
937:      Q NCswitch,NCcut,STT,SA,SB,SC) 
938:         USE KEY888:         USE KEY
939:         implicit NONE889:         implicit NONE
940:         integer I,N,J,AN,NATOMS,ATOMTYPES(NATOMS),MAXATOMTYPES890:         integer I, N, J, AN, NATOMS
941: 891: 
942: 892: 
943:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),893:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
944:      Q  grad(3*NATOMS), energy,STT(MAXATOMTYPES**2),ET894:      Q  grad(3*NATOMS), energy,STT,SA,SB,SC,ET
945:        DOUBLE PRECISION SA, SB, SC,SAT,SBT,SCT,STTT 
946:        DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES) 
947:        DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES) 
948:        DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES) 
949: 895: 
950: 896:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj
951:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj,POS 
952:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 897:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
953:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,MAXSEPSYS898:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,MAXSEPSYS
954:         LOGICAL NNCINC899:         LOGICAL NNCINC
955:         DIMENSION NNCINC(NATOMS*MAXSEP)900:         DIMENSION NNCINC(NATOMS*MAXSEP)
956: !$    INTEGER  OMP_GET_NUM_THREADS,901: !$    INTEGER  OMP_GET_NUM_THREADS,
957: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS902: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
958: 903: 
959:         INTEGER NNEXL1(NEXCLUSIONS)904:         INTEGER NNEXL1(NEXCLUSIONS)
960:         INTEGER NNEXL2(NEXCLUSIONS)905:         INTEGER NNEXL2(NEXCLUSIONS)
961:         DOUBLE PRECISION dx,dy,dz906:         DOUBLE PRECISION dx,dy,dz
974:         integer MaxGridX,MaxGridY,MaxGridZ919:         integer MaxGridX,MaxGridY,MaxGridZ
975:         double precision gridsize,RD1920:         double precision gridsize,RD1
976:         double precision minX,minY,minZ,maxX,maxY,maxZ921:         double precision minX,minY,minZ,maxX,maxY,maxZ
977:         integer Xgrid,Ygrid,Zgrid,TID922:         integer Xgrid,Ygrid,Zgrid,TID
978:         Rdiff=NCcut-NCswitch923:         Rdiff=NCcut-NCswitch
979:         alpha=12924:         alpha=12
980:         GRIDSIZE=NCcut*1.01925:         GRIDSIZE=NCcut*1.01
981:         Rcut2=NCcut**2926:         Rcut2=NCcut**2
982:         Rswitch2=NCswitch**2927:         Rswitch2=NCswitch**2
983: 928: 
984: !        SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 929:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
985: !     Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )930:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
986: 931: 
987: !        SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)932:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)
988: 933: 
989: !        SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)934:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)
990: ! grid the system935: ! grid the system
991:         minX=10000000936:         minX=10000000
992:         minY=10000000937:         minY=10000000
993:         minZ=10000000938:         minZ=10000000
994:         939:         
995:         maxX=-10000000940:         maxX=-10000000
996:         maxY=-10000000941:         maxY=-10000000
997:         maxZ=-10000000942:         maxZ=-10000000
998: 943: 
999:         do i=1,NATOMS944:         do i=1,NATOMS
1052:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i997:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1053:         enddo998:         enddo
1054: 999: 
1055: 1000: 
1056:            tempN=01001:            tempN=0
1057: 1002: 
1058: ! add a second loop that goes over charged atoms1003: ! add a second loop that goes over charged atoms
1059: 1004: 
1060: !$OMP PARALLEL1005: !$OMP PARALLEL
1061: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,1006: !$OMP& PRIVATE(i,ii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,
1062: !$OMP& Xgrid,Ygrid,Zgrid,TID,DINDEX,POS,SAT,SBT,SCT,STTT)  1007: !$OMP& Xgrid,Ygrid,Zgrid,TID,DINDEX)  
1063: !$OMP& REDUCTION(+:ENERGY,grad)1008: !$OMP& REDUCTION(+:ENERGY,grad)
1064:         TID=01009:         TID=0
1065:         NTHREADS=1;1010:         NTHREADS=1;
1066: !$      TID = OMP_GET_THREAD_NUM()1011: !$      TID = OMP_GET_THREAD_NUM()
1067: !$      NTHREADS = OMP_GET_NUM_THREADS()1012: !$      NTHREADS = OMP_GET_NUM_THREADS()
1068:         do i=1+TID,NATOMS,NTHREADS1013:         do i=1+TID,NATOMS,NTHREADS
1069: 1014: 
1070:           Xgrid=int((X(I)-minX)/gridsize)+11015:           Xgrid=int((X(I)-minX)/gridsize)+1
1071:           Ygrid=int((Y(I)-minY)/gridsize)+11016:           Ygrid=int((Y(I)-minY)/gridsize)+1
1072:           Zgrid=int((Z(I)-minZ)/gridsize)+11017:           Zgrid=int((Z(I)-minZ)/gridsize)+1
1073:           do ii=XGRID-1,XGRID+11018:           do ii=XGRID-1,XGRID+1
1074:            do jj=YGRID-1,YGRID+11019:            do jj=YGRID-1,YGRID+1
1075:             do kk=ZGRID-1,ZGRID+11020:             do kk=ZGRID-1,ZGRID+1
1076:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1021:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1077:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1022:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1078:               do jjj=1,Ngrid(ii,jj,kk)1023:               do jjj=1,Ngrid(ii,jj,kk)
1079:                C2=grid(ii,jj,kk,jjj)1024:                C2=grid(ii,jj,kk,jjj)
1080:                DINDEX=C2-I1025:                DINDEX=C2-I
1081:                if(DINDEX .GT. 0)THEN1026:                if(DINDEX .GT. 0)THEN
1082:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX)  ) )then1027:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX)  ) )then
 1028:                 if(I .EQ. 5178 .AND. C2 .EQ. 5212)THEN
 1029:                 endif
1083:                  dx = X(I) - X(C2)1030:                  dx = X(I) - X(C2)
1084:                  dy = Y(I) - Y(C2)1031:                  dy = Y(I) - Y(C2)
1085:                  dz = Z(I) - Z(C2)1032:                  dz = Z(I) - Z(C2)
1086:                  r2 = dx**2 + dy**2 + dz**21033:                  r2 = dx**2 + dy**2 + dz**2
1087:                  if(r2 .le. Rcut2)then1034:                  if(r2 .le. Rcut2)then
1088:                   POS=(ATOMTYPES(I)-1)*MAXATOMTYPES+ATOMTYPES(C2) 
1089:                   STTT=STT(POS) 
1090:                   SCT=SC(POS) 
1091:                   rm2 = 1/r21035:                   rm2 = 1/r2
1092:                   rm14 = rm2**71036:                   rm14 = rm2**7
1093:                   energy=energy+STTT*rm2**6+SCT1037:                   energy=energy+STT*rm2**6+SC
1094:                   f_over_r=-STTT*12.0*rm141038:                   f_over_r=-STT*12.0*rm14
1095:                   RD1=sqrt(r2)-NCswitch1039:                   RD1=sqrt(r2)-NCswitch
1096:                   if(r2 .gt. Rswitch2)then1040:                   if(r2 .gt. Rswitch2)then
1097:                    SAT=SA(POS)1041:                    f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
1098:                    SBT=SB(POS)1042:                    energy=energy+SA/3.0*RD1**3+SB/4.0*RD1**4
1099:                    f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2) 
1100:                    energy=energy+SAT/3.0*RD1**3+SBT/4.0*RD1**4 
1101:                   elseif(r2 .lt. Rswitch2)then1043:                   elseif(r2 .lt. Rswitch2)then
1102:                  ! normal repulsive term1044:                  ! normal repulsive term
1103:                   else1045:                   else
1104:                  ! things should have fallen in one of the previous two...1046:                  ! things should have fallen in one of the previous two...
1105:                    write(*,*) 'something went wrong with switching function'1047:                    write(*,*) 'something went wrong with switching function'
1106:                   endif1048:                   endif
1107:                   grad(I*3-2) = grad(I*3-2) + f_over_r * dx1049:                   grad(I*3-2) = grad(I*3-2) + f_over_r * dx
1108:                   grad(I*3-1) = grad(I*3-1) + f_over_r * dy1050:                   grad(I*3-1) = grad(I*3-1) + f_over_r * dy
1109:                   grad(I*3)   = grad(I*3)   + f_over_r * dz1051:                   grad(I*3)   = grad(I*3)   + f_over_r * dz
1110:                   grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx1052:                   grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
1122: 1064: 
1123: !$OMP DO1065: !$OMP DO
1124:        DO I=1,NEXCLUSIONS1066:        DO I=1,NEXCLUSIONS
1125:          C1=NNEXL1(I) 1067:          C1=NNEXL1(I) 
1126:          C2=NNEXL2(I)1068:          C2=NNEXL2(I)
1127:          dx = X(C1) - X(C2)1069:          dx = X(C1) - X(C2)
1128:          dy = Y(C1) - Y(C2)1070:          dy = Y(C1) - Y(C2)
1129:          dz = Z(C1) - Z(C2)1071:          dz = Z(C1) - Z(C2)
1130:          r2 = dx**2 + dy**2 + dz**21072:          r2 = dx**2 + dy**2 + dz**2
1131:          if(r2 .le. Rcut2)then1073:          if(r2 .le. Rcut2)then
1132:            POS=(ATOMTYPES(C1)-1)*MAXATOMTYPES+ATOMTYPES(C2)1074:           rm2 = 1/r2
1133:            STTT=STT(POS)1075:           rm14 = rm2**7
1134:            SCT=SC(POS)1076:            energy=energy-STT*rm2**6-SC
1135:  1077:            ET=-STT*rm2**6-SC
1136:            rm2 = 1/r21078:            f_over_r=-STT*12.0*rm14
1137:            rm14 = rm2**71079:            RD1=sqrt(r2)-NCswitch
1138:             energy=energy-STTT*rm2**6-SCT 
1139:             f_over_r=-STTT*12.0*rm14 
1140:             RD1=sqrt(r2)-NCswitch 
1141:           if(r2 .gt. Rswitch2)then1080:           if(r2 .gt. Rswitch2)then
1142:            SAT=SA(POS)1081:            f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
1143:            SBT=SB(POS)1082:            energy=energy-SA/3.0*RD1**3-SB/4.0*RD1**4
1144:            f_over_r=f_over_r+(SAT*RD1**2+SBT*RD1**3)*sqrt(rm2)1083:            ET=ET-SA/3.0*RD1**3-SB/4.0*RD1**4
1145:            energy=energy-SAT/3.0*RD1**3-SBT/4.0*RD1**4 
1146: !          elseif(r2 .lt. Rswitch2)then1084: !          elseif(r2 .lt. Rswitch2)then
1147:          ! normal repulsive term1085:          ! normal repulsive term
1148: !          else1086: !          else
1149:          ! things should have fallen in one of the previous two...1087:          ! things should have fallen in one of the previous two...
1150: !           write(*,*) 'something went wrong with switching function'1088: !           write(*,*) 'something went wrong with switching function'
1151:           endif1089:           endif
1152:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx1090:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx
1153:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy1091:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy
1154:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz1092:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz
1155:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx1093:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0