hdiff output

r30681/SBM.f 2016-07-06 15:38:22.513481811 +0100 r30680/SBM.f 2016-07-06 15:38:22.841486231 +0100
  2:       USE KEY  2:       USE KEY
  3:       implicit NONE  3:       implicit NONE
  4:       INTEGER NATOMS,i,j  4:       INTEGER NATOMS,i,j
  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)
  6:       DOUBLE PRECISION ENERGY,STT  6:       DOUBLE PRECISION ENERGY,STT
  7:   7: 
  8:       LOGICAL :: CALLED=.FALSE.  8:       LOGICAL :: CALLED=.FALSE.
  9:       LOGICAL GTEST, STEST  9:       LOGICAL GTEST, STEST
 10:         integer NSBMMAX 10:         integer NSBMMAX
 11:         parameter(NSBMMAX=25000) 11:         parameter(NSBMMAX=25000)
 12:       INTEGER MAXSEP,MAXSEPSYS 12:       INTEGER MAXSEP,MAXSEP1
 13:       PARAMETER (MAXSEP=20) 13:       PARAMETER (MAXSEP1=30)
 14:       LOGICAL NNCINC 14:       LOGICAL NNCINC(NSBMMAX,MAXSEP1)
 15:       DIMENSION NNCINC(NSBMMAX*MAXSEP) 
 16:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 15:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX),
 17:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF(NSBMMAX*5*6), 16:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF1(NSBMMAX*5),
 18:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 17:      Q CONCOEF2(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
 19:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2), 18:       INTEGER  Ib1(NSBMMAX*2), Ib2(NSBMMAX*2), IT(NSBMMAX*2), JT(NSBMMAX*2),
 20:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3), 19:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3),
 21:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5), 20:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5),
 22:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX) 21:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX)
 23:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NEXCLUDE,NEXCLUDEELEC 22:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NEXCLUDE,NEXCLUDEELEC
 24:       PARAMETER (MAXEXCLUSIONS=20) 23:       PARAMETER (MAXEXCLUSIONS=20)
 25:       PARAMETER (MAXEXCLUSIONSELEC=5) 24:       PARAMETER (MAXEXCLUSIONSELEC=5)
 26:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS), 25:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS),
 27:      Q NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC) 26:      Q NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC)
 28:  27: 
 29:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3) 28:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3)
 30:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX) 29:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX)
 31:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 30:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 32:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 31:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
 33:      Q SBMPRK,SBMPRX,CONCOEF 32:      Q SBMPRK,SBMPRX,CONCOEF1,CONCOEF2
 34:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC, 33:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC,
 35:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1, 34:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,
 36:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC 35:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC
 37:         common /logical/ NNCINC 36:         common /logical/ NNCINC
  37:         MAXSEP=MAXSEP1
 38:  38: 
 39:         if(NATOMS.gt. NSBMMAX)then 39:         if(NATOMS.gt. NSBMMAX)then
 40:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 40:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 41:         STOP 41:         STOP
 42:         endif 42:         endif
 43:  43: 
 44:        if(.NOT.CALLED)then 44:        if(.NOT.CALLED)then
 45:         NEXCLUDE=MAXEXCLUSIONS 45:         NEXCLUDE=MAXEXCLUSIONS
 46:         NEXCLUDEELEC=MAXEXCLUSIONS 46:         NEXCLUDEELEC=MAXEXCLUSIONS
 47:         MAXSEPSYS=0 
 48:         write(*,*) 47:         write(*,*)
 49:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 48:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 50:         write(*,*)  '   For more information on SMOG models, see:' 49:         write(*,*)  '   For more information on SMOG models, see:'
 51:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 50:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 52:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 51:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 53:         write(*,*) 52:         write(*,*)
 54:          53:         
 55:  54: 
 56:        call SBMinit(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP, 55:        call SBMinit(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,
 57:      Q PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE, 56:      Q PK,PHITYPE,PHISBM,IC,JC,CONCOEF1,CONCOEF2,NNEXL1,NNEXL2,NEXCLUDE,
 58:      Q NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC, 57:      Q NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,NBA,NTA,NPA,NC,
 59:      Q SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,NSBMMAX,CONTACTTYPE, 58:      Q SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,NSBMMAX,CONTACTTYPE,
 60:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 59:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 61:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or. 60:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or.
 62:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then 61:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then
 63:         write(*,*) 'increase array size' 62:         write(*,*) 'increase array size'
 64:         write(*,*) NBA, NTA,NPA,NC  63:         write(*,*) NBA, NTA,NPA,NC 
 65:         stop 64:         stop
 66:         endif  65:         endif 
 67:         CALLED=.TRUE. 66:         CALLED=.TRUE.
 68:         endIF 67:         endIF
 69:  68: 
 70:  69: 
 71: ! call the energy routine 70: ! call the energy routine
 72:       call calc_energy_SBM(qo,natoms,NSBMMAX,GRAD,energy,Ib1,Ib2,Rb,bK,IT,JT,KT, 71:       call calc_energy_SBM(qo,natoms,GRAD,energy,Ib1,Ib2,Rb,bK,IT,JT,KT,
 73:      Q ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF,  72:      Q ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF1,CONCOEF2, 
 74:      Q NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC, 73:      Q NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,
 75:      Q MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,STT,NCcut, 74:      Q MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,STT,NCcut,
 76:      Q CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK, 75:      Q CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,
 77:      Q SBMPRX) 76:      Q SBMPRX)
 78:       IF (STEST) THEN 77:       IF (STEST) THEN
 79:          PRINT '(A)','ERROR - second derivatives not available' 78:          PRINT '(A)','ERROR - second derivatives not available'
 80:          STOP 79:          STOP
 81:       ENDIF 80:       ENDIF
 82:       return 81:       return
 83:       end 82:       end
 84:  83: 
 85:  84: 
 86: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 85: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 87: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 86: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 88: !* startt then the velocities are assigned, otherwise, they are read   * 87: !* startt then the velocities are assigned, otherwise, they are read   *
 89: !* by selecting 2, or generated by selecting 3                         * 88: !* by selecting 2, or generated by selecting 3                         *
 90: !*********************************************************************** 89: !***********************************************************************
 91:  90: 
 92:       subroutine SBMINIT(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk, 91:       subroutine SBMINIT(NATOMS,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,
 93:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF,  92:      Q IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC, JC, CONCOEF1,CONCOEF2, 
 94:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNEXL1ELEC,NNEXL2ELEC,MAXEXCLUSIONSELEC, 93:      Q NNEXL1,NNEXL2,MAXEXCLUSIONS,NNEXL1ELEC,NNEXL2ELEC,MAXEXCLUSIONSELEC,
 95:      Q NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON, 94:      Q NNCINC,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,
 96:      Q NCswitch,NCcut,STT,NSBMMAX,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch, 95:      Q STT,NSBMMAX,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,
 97:      Q DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 96:      Q SBMPRI,SBMPRK,SBMPRX)
 98:       USE KEY 97:       USE KEY
 99:       USE COMMONS, only: ATMASS 98:       USE COMMONS, only: ATMASS
100:       implicit NONE 99:       implicit NONE
101: 100: 
102:         integer i,j,K,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,101:         integer i,j,MaxCon,NATOMS,storage, dummy,ANr, Ib22,Ib21,IT1,JT1,
103:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 102:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 
104:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,NSBMMAX,103:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,NSBMMAX,
105:      Q MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,MAXSEPSYS104:      Q MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,MAXSEPTEMP
106: 105: 
107:       DOUBLE PRECISION :: Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 106:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX), 
108:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX),107:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX),
109:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS)108:      Q Sigma, EpsC, NCswitch,NCcut,STT,SBMCHARGE(NATOMS),
110: !      DOUBLE PRECISION CONCOEF109:      Q CONCOEF1(5*NSBMMAX),CONCOEF2(5*NSBMMAX) 
111:       DOUBLE PRECISION CONCOEF (NC*6) 
112:       INTEGER  Ib1(2*NSBMMAX), Ib2(2*NSBMMAX), IT(2*NSBMMAX), JT(2*NSBMMAX), 110:       INTEGER  Ib1(2*NSBMMAX), Ib2(2*NSBMMAX), IT(2*NSBMMAX), JT(2*NSBMMAX), 
113:      Q KT(2*NSBMMAX),IP(3*NSBMMAX),JP(3*NSBMMAX),KP(3*NSBMMAX),111:      Q KT(2*NSBMMAX),IP(3*NSBMMAX),JP(3*NSBMMAX),KP(3*NSBMMAX),
114:      Q LP(3*NSBMMAX),IC(NSBMMAX*5),JC(NSBMMAX*5), 112:      Q LP(3*NSBMMAX),IC(NSBMMAX*5),JC(NSBMMAX*5), 
115:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE113:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE
116:       INTEGER PHITYPE(3*NSBMMAX),M114:       INTEGER PHITYPE(3*NSBMMAX),M
117:       integer AA,BB,ANTEMP115:       integer AA,BB,ANTEMP
118:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)116:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)
119:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)117:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)
120:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)118:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)
121:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)119:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)
122:       DOUBLE PRECISION dx,dy,dz120:       DOUBLE PRECISION dx,dy,dz
123:       double precision PI121:       double precision PI
124:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut122:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
125:       integer CONTACTTYPE(NSBMMAX*5),TMPARG123:       integer CONTACTTYPE(NSBMMAX*5)
126:       character TMPCHAR124:       character TMPCHAR
127:       integer TMPINT,NUMCHARGES,nexc,I1,I2125:       integer TMPINT,NUMCHARGES,nexc,I1,I2
128:       double precision TMPREAL,concentration126:       double precision TMPREAL,concentration
129:       LOGICAL NNCINC(NATOMS*MAXSEP)127:       LOGICAL NNCINC(NATOMS,MAXSEP)
130: 128: 
131:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)129:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)
132:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)130:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
133:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS131:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS
134:       DIMENSION EXCLUSIONS(NATOMS,MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)132:       DIMENSION EXCLUSIONS(NATOMS,MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)
135: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 133: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
136: !$OMP PARALLEL134: !$OMP PARALLEL
137: !$    NTHREADS = OMP_GET_NUM_THREADS()135: !$    NTHREADS = OMP_GET_NUM_THREADS()
138: !$      TID = OMP_GET_THREAD_NUM()136: !$      TID = OMP_GET_THREAD_NUM()
139: !$    if(TID .eq. 0)then 137: !$    if(TID .eq. 0)then 
140: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS138: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS
141: !$    endif139: !$    endif
142: !$OMP END PARALLEL140: !$OMP END PARALLEL
143: 141: 
 142: 
144:       DO I=1,NATOMS143:       DO I=1,NATOMS
145:         DO J=1,MAXSEP144:         DO J=1,MAXSEP
146:           NNCINC(I*MAXSEP-J+1)=.TRUE.145:           NNCINC(I,J)=.TRUE.
147:         enddo146:         enddo
148:       enddo147:       enddo
149: 148: 
150:       pi = 3.14159265358979323846264338327950288419716939937510149:       pi = 3.14159265358979323846264338327950288419716939937510
151:       MaxCon=NATOMS*5150:       MaxCon=NATOMS*5
 151:       MAXSEPTEMP=0
152:       do i=1,NATOMS152:       do i=1,NATOMS
153:         NUMOFEXCLUSIONS(I)=0153:         NUMOFEXCLUSIONS(I)=0
154:       enddo154:       enddo
155: 155: 
156: ! These lines read in the parameters.156: ! These lines read in the parameters.
157:         open(30, file='SBM.INP', status='old', access='sequential')157:         open(30, file='SBM.INP', status='old', access='sequential')
158:         write(*,*) 158:         write(*,*) 
159:         write(*,*) 'Reading the forcefield from SBM.INP'159:         write(*,*) 'Reading the forcefield from SBM.INP'
160:         write(*,*) 160:         write(*,*) 
161:         read(30,*)161:         read(30,*)
187:         read(30,*) NC187:         read(30,*) NC
188:         write(*,*) NC, 'contacts'        188:         write(*,*) NC, 'contacts'        
189:           if(NC .gt. MaxCon)then189:           if(NC .gt. MaxCon)then
190:              write(*,*) 'too many contacts'190:              write(*,*) 'too many contacts'
191:              STOP191:              STOP
192:           endif192:           endif
193: 193: 
194:         do i=1, NC194:         do i=1, NC
195:           read(30, *) IC(i), JC(i), CONTACTTYPE(I),Sigma, EpsC195:           read(30, *) IC(i), JC(i), CONTACTTYPE(I),Sigma, EpsC
196:           !! Add more contact types196:           !! Add more contact types
197:           TMPARG=0197:           CALL INCLUDEEXCLUSIONS(NATOMS,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)
198:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,TMPARG) 
199:           if(CONTACTTYPE(I) .eq. 1)then198:           if(CONTACTTYPE(I) .eq. 1)then
200:             CONCOEF(I*6-5)=EPSC*SIGMA**12199:             ConCoef1(I)=EPSC*SIGMA**12-REPS*RSIG**12
201:             CONCOEF(I*6-4)=2*EPSC*SIGMA**6200:             ConCoef2(I)=2*EPSC*SIGMA**6
202:           elseif(CONTACTTYPE(I) .eq. 2)then201:           elseif(CONTACTTYPE(I) .eq. 2)then
203:             CONCOEF(I*6-5)=5*EPSC*SIGMA**12202:             ConCoef1(I)=5*EPSC*SIGMA**12-REPS*RSIG**12
204:             CONCOEF(I*6-4)=6*EPSC*SIGMA**10203:             ConCoef2(I)=6*EPSC*SIGMA**10
205:           else204:           else
206:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.'205:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.'
207:           endif206:           endif
208:         end do207:         end do
209: ! make a function to check and add208: ! make a function to check and add
210: 209: 
211:           read(30,*)210:           read(30,*)
212:           read(30,*) nBA211:           read(30,*) nBA
213:           write(*,*) nBA, 'bonds'212:           write(*,*) nBA, 'bonds'
214:         do i=1, nBA213:         do i=1, nBA
215:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)214:           read(30,*) Ib1(i), Ib2(i),bondtype,Rb(i), bK(i)
216:           !!! ADD BONDS OF TYPE 6215:           !!! ADD BONDS OF TYPE 6
217:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP)THEN216:           IF(max(IB1(I),IB2(I))-min(IB1(I),IB2(I)) .gt. MAXSEP)THEN
218:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS LESS THAN MAXSEP'217:             WRITE(*,*) 'WARNING: DIFFERENCE IN BONDED ATOM INDICES IS LESS THAN MAXSEP'
219:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'218:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
220:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'219:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
221:           ENDIF220:           ENDIF
222:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN221:           IF(bondtype .eq. 1 .or. bondtype .eq. 6)THEN
223:             CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          222:             CALL INCLUDEEXCLUSIONS(NATOMS,IB1(I),IB2(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
224:           ENDIF223:           ENDIF
225:         end do224:         end do
226: 225: 
227:           read(30,*)226:           read(30,*)
228:           read(30,*) nTA227:           read(30,*) nTA
229:           write(*,*) nTA, 'angles'228:           write(*,*) nTA, 'angles'
230:         do i=1, nTA229:         do i=1, nTA
231:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)230:           read(30,*) IT(i), JT(i), KT(i), ANTC(i), Tk(i)
232:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          231:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),JT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
233:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          232:           CALL INCLUDEEXCLUSIONS(NATOMS,IT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
234:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          233:           CALL INCLUDEEXCLUSIONS(NATOMS,JT(I),KT(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
235:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.234:           IF(max(IT(I),JT(I))-min(IT(I),JT(I)) .gt. MAXSEP .OR.
236:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.235:      Q       max(IT(I),KT(I))-min(IT(I),KT(I)) .gt. MAXSEP .OR.
237:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN236:      Q       max(KT(I),JT(I))-min(KT(I),JT(I)) .gt. MAXSEP )THEN
238:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS LESS THAN MAXSEP'237:             WRITE(*,*) 'WARNING: DIFFERENCE IN BOND ANGLE ATOM INDICES IS LESS THAN MAXSEP'
239:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'238:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
240:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'239:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
241:           ENDIF240:           ENDIF
242:         ENDDO241:         ENDDO
243: 242: 
244:           read(30,*) 243:           read(30,*) 
245:           read(30,*) nPA244:           read(30,*) nPA
246:           write(*,*) nPA, 'dihedrals'245:           write(*,*) nPA, 'dihedrals'
247: 246: 
248: ! this reads in the dihedral angles and calculates the cosines and sines247: ! this reads in the dihedral angles and calculates the cosines and sines
249: ! in order to make the force and energy calculations easier, later.248: ! in order to make the force and energy calculations easier, later.
250:         do i=1, npA249:         do i=1, npA
251:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)250:           read(30,*) IP(i),JP(i),KP(i),LP(i),PHITYPE(i),PHISBM(i),PK(i)
252:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IP(I),JP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          251:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),JP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
253:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          252:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
254:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,IP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          253:           CALL INCLUDEEXCLUSIONS(NATOMS,IP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
255:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          254:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),KP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
256:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          255:           CALL INCLUDEEXCLUSIONS(NATOMS,JP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
257:           CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          256:           CALL INCLUDEEXCLUSIONS(NATOMS,KP(I),LP(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
258: 257: 
259:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.258:           IF(max(IP(I),JP(I))-min(IP(I),JP(I)) .gt. MAXSEP .OR.
260:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.259:      Q       max(IP(I),KP(I))-min(IP(I),KP(I)) .gt. MAXSEP .OR.
261:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.260:      Q       max(IP(I),LP(I))-min(IP(I),LP(I)) .gt. MAXSEP .OR.
262:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.261:      Q       max(JP(I),KP(I))-min(JP(I),KP(I)) .gt. MAXSEP .OR.
263:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.262:      Q       max(JP(I),LP(I))-min(JP(I),LP(I)) .gt. MAXSEP .OR.
264:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN263:      Q       max(KP(I),LP(I))-min(KP(I),LP(I)) .gt. MAXSEP )THEN
265:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS LESS THAN MAXSEP'264:             WRITE(*,*) 'WARNING: DIFFERENCE IN DIHEDRAL ANGLE ATOM INDICES IS LESS THAN MAXSEP'
266:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'265:             WRITE(*,*) '         THIS CAN LEAD TO POOR PERFORMANCE.'
267:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'266:             WRITE(*,*) '         CONSIDER RECOMPILING WITH A LARGER VALUE OF MAXSEP.'
268:           ENDIF267:           ENDIF
269:         ENDDO 268:         ENDDO 
270: 269: 
271:         read(30,*)270:         read(30,*)
272:         read(30,*) nexc271:         read(30,*) nexc
273:         write(*,*) nexc, 'exculusions'272:         write(*,*) nexc, 'exculusions'
274: 273: 
275:       do i=1, nexc274:       do i=1, nexc
276:         read(30,*) I1, I2275:         read(30,*) I1, I2
277:         CALL INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,I1,I2,EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS)          276:         CALL INCLUDEEXCLUSIONS(NATOMS,I1,I2,EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)          
278:       enddo277:       enddo
279: 278: 
280: ! now that we know what to exclude, let's make a static list of pairs to279: ! now that we know what to exclude, let's make a static list of pairs to
281: ! exclude280: ! exclude
282: ! At this point, we will change the meaning of maxexclusions, which is281: ! At this point, we will change the meaning of maxexclusions, which is
283: ! actually NEXCLUDE282: ! actually NEXCLUDE
284:         NT=MAXEXCLUSIONS*NATOMS      283:         NT=MAXEXCLUSIONS*NATOMS      
285:         NTE=MAXEXCLUSIONSELEC*NATOMS      284:         NTE=MAXEXCLUSIONSELEC*NATOMS      
286:         MAXEXCLUSIONS=0  285:         MAXEXCLUSIONS=0  
287:         MAXEXCLUSIONSELEC=0  286:         MAXEXCLUSIONSELEC=0  
309: 308: 
310:         do I=1,NATOMS309:         do I=1,NATOMS
311:           TARR(I)=0310:           TARR(I)=0
312:         enddo311:         enddo
313: 312: 
314: ! read in position restraints313: ! read in position restraints
315:        read(30,*) 314:        read(30,*) 
316:        read(30,*) SBMPRN315:        read(30,*) SBMPRN
317:        write(*,*) SBMPRN, 'position restraints'316:        write(*,*) SBMPRN, 'position restraints'
318:        do I=1,SBMPRN317:        do I=1,SBMPRN
319:             J=(I-1)*6318:             read(30,*) SBMPRI(I),SBMPRK(I,1),SBMPRK(I,2),SBMPRK(I,3),
320:             K=J/2319:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6),
321:             read(30,*) SBMPRI(I),SBMPRK(J+1),SBMPRK(J+2),SBMPRK(J+3),320:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)
322:      Q         SBMPRK(J+4),SBMPRK(J+5),SBMPRK(J+6), 
323:      Q         SBMPRX(K+1),SBMPRX(K+2),SBMPRX(K+3) 
324:             if(TARR(SBMPRI(I)) .ne. 0)then321:             if(TARR(SBMPRI(I)) .ne. 0)then
325:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)322:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
326:                 STOP323:                 STOP
327:             endif324:             endif
328:             TARR(SBMPRI(I))=1325:             TARR(SBMPRI(I))=1
329:         if(SBMPRK(J+4) .ne. 0 .or. SBMPRK(J+5) .ne. 0 .or. SBMPRK(J+6) .ne.0)then  326:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  
330:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'327:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'
331:         endif328:         endif
332:        enddo 329:        enddo 
 330:         MAXSEP=MAXSEPTEMP
333:        close(30)331:        close(30)
334:       END 332:       END 
335: 333: 
336: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^334: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMinit^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
337: 335: 
338: 336: 
339: !**************************************************************************337: !**************************************************************************
340: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in338: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in
341: ! the model339: ! the model
342: !**************************************************************************340: !**************************************************************************
343: 341: 
344:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,NSBMMAX,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,342:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,
345:      Q NNCINC,MAXSEP,MAXSEPSYS)343:      Q NNCINC,MAXSEP,MAXSEPTEMP)
346:       IMPLICIT NONE344:       IMPLICIT NONE
347:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,345:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,MAXSEPTEMP
348:      Q MAXSEPSYS,NSBMMAX 
349:       DIMENSION EXCLIST(NATOMS,MAXEXCLUDE),NEXCLUDE(NATOMS)346:       DIMENSION EXCLIST(NATOMS,MAXEXCLUDE),NEXCLUDE(NATOMS)
350:       LOGICAL  NNCINC347:       LOGICAL NNCINC(NATOMS,MAXSEP)
351:       DIMENSION NNCINC(NATOMS*MAXSEP) 
352:       INTEGER I,N,M,AA,BB,AB348:       INTEGER I,N,M,AA,BB,AB
353:       LOGICAL INCLUDED349:       LOGICAL INCLUDED
354: 350: 
355: ! If the atoms are within MAXSEP of each other, then keep track via351: ! If the atoms are within MAXSEP of each other, then keep track via
356: ! NNCINC352: ! NNCINC
357:       AA=min(ATOM1,ATOM2)353:       AA=min(ATOM1,ATOM2)
358:       BB=max(ATOM1,ATOM2)354:       BB=max(ATOM1,ATOM2)
359:       AB=BB-AA355:       AB=BB-AA
360:       IF(AB .le. MAXSEP)THEN356:       IF(AB .le. MAXSEP)THEN
361:         NNCINC((AA-1)*MAXSEP+AB)=.FALSE.357:         NNCINC(AA, AB)=.FALSE.
362:         IF(AB .gt. MAXSEPSYS)THEN358:         IF(AB .gt. MAXSEPTEMP)THEN
363:           MAXSEPSYS=AB359:           MAXSEPTEMP=AB
364:         ENDIF360:         ENDIF
365:       ELSE361:       ELSE
366:         INCLUDED = .FALSE.362:         INCLUDED = .FALSE.
367:         N=MAX(ATOM1,ATOM2)363:         N=MAX(ATOM1,ATOM2)
368:         M=MIN(ATOM1,ATOM2)364:         M=MIN(ATOM1,ATOM2)
369:         DO I =1,NEXCLUDE(M)365:         DO I =1,NEXCLUDE(M)
370:           IF(EXCLIST(M,I) .eq. N)THEN366:           IF(EXCLIST(M,I) .eq. N)THEN
371:             INCLUDED = .TRUE.367:             INCLUDED = .TRUE.
372:           ENDIF368:           ENDIF
373:         ENDDO 369:         ENDDO 
380:           EXCLIST(M,NEXCLUDE(M))=N376:           EXCLIST(M,NEXCLUDE(M))=N
381:         ENDIF377:         ENDIF
382:       ENDIF378:       ENDIF
383:       END379:       END
384: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^380: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
385: 381: 
386: 382: 
387: C383: C
388: C Calculate the Forces and energies384: C Calculate the Forces and energies
389: C385: C
390:       subroutine CALC_ENERGY_SBM(qo,natoms,NSBMMAX,GRAD,energy,Ib1,Ib2,386:       subroutine CALC_ENERGY_SBM(qo,natoms,GRAD,energy,Ib1,Ib2,Rb,bK,IT,
391:      Q Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,387:      Q JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF1,
392:      Q CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,NNEXL2ELEC,388:      Q CONCOEF2,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,NNEXL2ELEC,
393:      Q NEXCLUSIONSELEC,NCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,389:      Q NEXCLUSIONSELEC,NCINC,MAXSEP,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,
394:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,390:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
395:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)391:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
396: 392: 
397:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,393:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,NEXCLUSIONSELEC
398:      Q NSBMMAX,NEXCLUSIONSELEC 
399: 394: 
400:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY395:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
401:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)396:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
402:       LOGICAL NNCINC397:       LOGICAL NNCINC(NATOMS,MAXSEP)
403:       DIMENSION NNCINC(NATOMS*MAXSEP) 
404: 398: 
405:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 399:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
406:      Q PHISBM(NPA),NCswitch,NCcut,STT, 400:      Q PHISBM(NPA),CONCOEF1(NC),CONCOEF2(NC),NCswitch,NCcut,STT, 
407:      Q SBMCHARGE(NATOMS)401:      Q SBMCHARGE(NATOMS)
408:         DOUBLE PRECISION  CONCOEF(NC*6) 
409:         INTEGER Ib1(NBA),Ib2(NBA),IT(NTA),JT(NTA),KT(NTA),IP(NPA), 402:         INTEGER Ib1(NBA),Ib2(NBA),IT(NTA),JT(NTA),KT(NTA),IP(NPA), 
410:      Q JP(NPA),KP(NPA),LP(NPA),IC(NC),JC(NC),SBMCHARGEON(NATOMS),403:      Q JP(NPA),KP(NPA),LP(NPA),IC(NC),JC(NC),SBMCHARGEON(NATOMS),
411:      Q PHITYPE(NPA),CONTACTTYPE(NC)404:      Q PHITYPE(NPA),CONTACTTYPE(NC)
412:         INTEGER NNEXL1(NEXCLUSIONS)405:         INTEGER NNEXL1(NEXCLUSIONS)
413:         INTEGER NNEXL2(NEXCLUSIONS)406:         INTEGER NNEXL2(NEXCLUSIONS)
414:         INTEGER NNEXL1ELEC(NEXCLUSIONSELEC)407:         INTEGER NNEXL1ELEC(NEXCLUSIONSELEC)
415:         INTEGER NNEXL2ELEC(NEXCLUSIONSELEC)408:         INTEGER NNEXL2ELEC(NEXCLUSIONSELEC)
416:       DOUBLE PRECISION dx,dy,dz409:       DOUBLE PRECISION dx,dy,dz
417:       INTEGER SBMPRN, SBMPRI(NATOMS)410:       INTEGER SBMPRN, SBMPRI(NATOMS)
418:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)411:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
419:       do i = 1, natoms412:       do i = 1, natoms
420:         j = (i-1)*3413:         j = (i-1)*3
421:         x(i) = qo(j+1)414:         x(i) = qo(j+1)
422:         y(i) = qo(j+2)415:         y(i) = qo(j+2)
423:         z(i) = qo(j+3)416:         z(i) = qo(j+3)
424:         grad(j+1) = 0.0417:         grad(j+1) = 0.0
425:         grad(j+2) = 0.0418:         grad(j+2) = 0.0
426:         grad(j+3) = 0.0419:         grad(j+3) = 0.0
427:       enddo420:       enddo
428: 421: 
429:       energy = 0.0422:       energy = 0.0
430:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)423:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
431:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)424:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
432:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,425:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
433:      Q PHITYPE,PHISBM,NPA)426:      Q PHITYPE,PHISBM,NPA)
434:       call SBMContacts(x,y,z,grad, energy,natoms,NSBMMAX,IC,JC,CONCOEF,427:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC,CONCOEF1,CONCOEF2,
435:      Q NC,CONTACTTYPE)428:      Q NC,CONTACTTYPE)
436:       call SBMNonContacts(x,y,z,grad,energy,natoms,NSBMMAX,NNEXL1,NNEXL2,429:       call SBMNonContacts(x,y,z,grad, energy, natoms,NNEXL1,NNEXL2,
437:      Q NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)430:      Q NEXCLUSIONS,NNCINC,MAXSEP,NCswitch,NCcut,STT)
438:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,431:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,
439:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)432:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
440:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)433:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
441: 434: 
442:       end435:       end
443: 436: 
444: 437: 
445: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<438: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
446: !* SBMBonds  computes the hookean force and energy between chosen atoms *439: !* SBMBonds  computes the hookean force and energy between chosen atoms *
447: !***********************************************************************440: !***********************************************************************
734: 727: 
735: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^728: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^
736: 729: 
737: 730: 
738: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<731: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
739: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *732: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *
740: !* 10-12 or 6-12 potential                                              *733: !* 10-12 or 6-12 potential                                              *
741: !***********************************************************************734: !***********************************************************************
742: 735: 
743:       subroutine SBMcontacts(x,y,z,grad,energy,736:       subroutine SBMcontacts(x,y,z,grad,energy,
744:      Q NATOMS,NSBMMAX,IC,JC,CONCOEF,NC,CONTACTTYPE)737:      Q NATOMS,IC,JC,CONCOEF1,CONCOEF2,NC,CONTACTTYPE)
745:       USE KEY738:       USE KEY
746:       implicit NONE739:       implicit NONE
747:       integer I, NATOMS,NC,NSBMMAX740:       integer I, NATOMS,NC
748: 741: 
749:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 742:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
750:      Q grad(3*NATOMS), energy743:      Q grad(3*NATOMS), energy
751:       DOUBLE PRECISION dx,dy,dz,CO1,CO2744:       DOUBLE PRECISION dx,dy,dz
752: 745: 
753:       integer C1, C2 746:       integer C1, C2 
754:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 747:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 
755: 748: 
756:         DOUBLE PRECISION CONCOEF(NC*6)749:         DOUBLE PRECISION CONCOEF1(NC), CONCOEF2(NC)
757:         INTEGER IC(NC), JC(NC),CONTACTTYPE(NC)750:         INTEGER IC(NC), JC(NC),CONTACTTYPE(NC)
758: 751: 
759: ! type 1 is 6-12752: ! type 1 is 6-12
760: ! type 2 is 12-12753: ! type 2 is 12-12
761: ! implement type 3 as gaussian754: ! implement type 3 as gaussian
762: ! implement type 4 as dual-basin gaussian755: ! implement type 4 as dual-basin gaussian
763: ! question: Should type be part of the interaction, or should they be organized756: ! question: Should type be part of the interaction, or should they be organized
764: ! into separate lists?757: ! into separate lists?
 758: 
765: ! type 1 is 6-12 interaction759: ! type 1 is 6-12 interaction
766: !$OMP PARALLEL PRIVATE(I,CO1,CO2,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R) REDUCTION(+:ENERGY,grad)760: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R) REDUCTION(+:ENERGY,grad)
767: !$OMP DO761: !$OMP DO
768:           do i=1, NC762:           do i=1, NC
769:             C1 = IC(i)763:             C1 = IC(i)
770:             C2 = JC(i)764:             C2 = JC(i)
771:             dx = X(C1) - X(C2)765:             dx = X(C1) - X(C2)
772:             dy = Y(C1) - Y(C2)766:             dy = Y(C1) - Y(C2)
773:             dz = Z(C1) - Z(C2)767:             dz = Z(C1) - Z(C2)
774:             r2 = dx**2 + dy**2 + dz**2768:             r2 = dx**2 + dy**2 + dz**2
775:             rm2 = 1.0/r2769:             rm2 = 1.0/r2
776: 770: 
777:             if(CONTACTTYPE(I) .eq. 1)then771:             if(CONTACTTYPE(I) .eq. 1)then
778:               ! type 1 is 6-12 interaction772:               ! type 1 is 6-12 interaction
779:               rm6 = rm2**3773:               rm6 = rm2**3
780:               CO1=CONCOEF(I*6-5)774:               ENERGY = ENERGY + rm6*(CONCOEF1(I)*rm6-CONCOEF2(I))
781:               CO2=CONCOEF(I*6-4)775:               f_over_r = -12.0*rm6*(CONCOEF1(I)*rm6-CONCOEF2(I)/2.0)*rm2
782:               ENERGY = ENERGY + rm6*(CO1*rm6-CO2) 
783:               f_over_r = -12.0*rm6*(CO1*rm6-CO2/2.0)*rm2 
784:             elseif(CONTACTTYPE(I) .eq. 2)then776:             elseif(CONTACTTYPE(I) .eq. 2)then
785:               ! type 2 is 10-12 interaction777:               ! type 2 is 10-12 interaction
786:               rm10 = rm2**5778:               rm10 = rm2**5
787:               CO1=CONCOEF(I*6-5)779:               ENERGY = ENERGY + rm10*(CONCOEF1(I)*rm2-CONCOEF2(I))
788:               CO2=CONCOEF(I*6-4)780:               f_over_r = -60.0*rm10*(CONCOEF1(I)*rm2/5.0-CONCOEF2(I)/6.0)*rm2
789:               ENERGY = ENERGY + rm10*(CO1*rm2-CO2) 
790:               f_over_r = -60.0*rm10*(CO1*rm2/5.0-CO2/6.0)*rm2 
791:             else781:             else
792:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)782:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I)
793:               STOP783:               STOP
794:             endif784:             endif
795: 785: 
796:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx786:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
797:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy787:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
798:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz788:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz
799:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx789:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
800:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy790:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
806: 796: 
807:       end797:       end
808: 798: 
809: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^799: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
810: 800: 
811: 801: 
812: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<802: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
813: !* SBMNonContacts computes the forces due to non native contacts       *803: !* SBMNonContacts computes the forces due to non native contacts       *
814: !**********************************************************************804: !**********************************************************************
815: 805: 
816:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,NSBMMAX,NNEXL1,NNEXL2,806:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,NNEXL1,NNEXL2,
817:      Q  NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,NCswitch,NCcut,STT)807:      Q  NEXCLUSIONS,NNCINC,MAXSEP,NCswitch,NCcut,STT)
818:         USE KEY808:         USE KEY
819:         implicit NONE809:         implicit NONE
820:         integer I, N, J, AN, NATOMS810:         integer I, N, J, AN, NATOMS
821: 811: 
822: 812: 
823:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),813:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
824:      Q  grad(3*NATOMS), energy,STT,SA,SB,SC814:      Q  grad(3*NATOMS), energy,STT,SA,SB,SC
825: 815: 
826:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj816:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj
827:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 817:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
828:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,NSBMMAX,MAXSEPSYS818:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP
829:         LOGICAL NNCINC819:         LOGICAL NNCINC(NATOMS,MAXSEP)
830:         DIMENSION NNCINC(NATOMS*MAXSEP) 
831: !$    INTEGER  OMP_GET_NUM_THREADS,820: !$    INTEGER  OMP_GET_NUM_THREADS,
832: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS821: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
833: 822: 
834:         INTEGER NNEXL1(NEXCLUSIONS)823:         INTEGER NNEXL1(NEXCLUSIONS)
835:         INTEGER NNEXL2(NEXCLUSIONS)824:         INTEGER NNEXL2(NEXCLUSIONS)
836:         DOUBLE PRECISION dx,dy,dz825:         DOUBLE PRECISION dx,dy,dz
837:         integer tempN, alpha826:         integer tempN, alpha
838:         double precision Rdiff,Vfunc,Ffunc827:         double precision Rdiff,Vfunc,Ffunc
839:         double precision Rcut2,Rswitch2828:         double precision Rcut2,Rswitch2
840:         ! Ngrid is the number of atoms in that grid point829:         ! Ngrid is the number of atoms in that grid point
845:         ! dimensions of grid834:         ! dimensions of grid
846:         parameter (maxgrid=30)835:         parameter (maxgrid=30)
847:         dimension Ngrid(maxgrid,maxgrid,maxgrid),836:         dimension Ngrid(maxgrid,maxgrid,maxgrid),
848:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)837:      Q  grid(maxgrid,maxgrid,maxgrid,maxpergrid)
849:         integer MaxGridX,MaxGridY,MaxGridZ838:         integer MaxGridX,MaxGridY,MaxGridZ
850:         double precision gridsize,RD1839:         double precision gridsize,RD1
851:         double precision minX,minY,minZ,maxX,maxY,maxZ840:         double precision minX,minY,minZ,maxX,maxY,maxZ
852:         integer Xgrid,Ygrid,Zgrid,TID841:         integer Xgrid,Ygrid,Zgrid,TID
853:         Rdiff=NCcut-NCswitch842:         Rdiff=NCcut-NCswitch
854:         alpha=12843:         alpha=12
 844: 
855:         GRIDSIZE=NCcut*1.01845:         GRIDSIZE=NCcut*1.01
856:         Rcut2=NCcut**2846:         Rcut2=NCcut**2
857:         Rswitch2=NCswitch**2847:         Rswitch2=NCswitch**2
858: 848: 
859:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 849:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
860:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )850:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
861: 851: 
862:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)852:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)
863: 853: 
864:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)854:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)
949:            do jj=YGRID-1,YGRID+1939:            do jj=YGRID-1,YGRID+1
950:             do kk=ZGRID-1,ZGRID+1940:             do kk=ZGRID-1,ZGRID+1
951:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.941:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
952:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then942:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
953:               do jjj=1,Ngrid(ii,jj,kk)943:               do jjj=1,Ngrid(ii,jj,kk)
954: 944: 
955:                C2=grid(ii,jj,kk,jjj)945:                C2=grid(ii,jj,kk,jjj)
956:                DINDEX=C2-I946:                DINDEX=C2-I
957:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNEXL(C1,C12))) )then947:            !if(C12 .gt. 0 .and. (C12 .gt. MAXSEP .or. (C12 .le. MAXSEP .and. NNEXL(C1,C12))) )then
958:            !if(DINDEX .gt. MAXSEP .or. (DINDEX .gt. 0 .and. NNEXL(I,DINDEX)))then948:            !if(DINDEX .gt. MAXSEP .or. (DINDEX .gt. 0 .and. NNEXL(I,DINDEX)))then
959:                if(DINDEX .gt. 0 .AND. ( DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX) ) ) )then949:                if(DINDEX .gt. 0 .AND. ( DINDEX .GT. MAXSEP .OR. (DINDEX .LE. MAXSEP .AND. NNCINC(I,DINDEX) ) ) )then
960:                 dx = X(I) - X(C2)950:                 dx = X(I) - X(C2)
961:                 dy = Y(I) - Y(C2)951:                 dy = Y(I) - Y(C2)
962:                 dz = Z(I) - Z(C2)952:                 dz = Z(I) - Z(C2)
963:                 r2 = dx**2 + dy**2 + dz**2953:                 r2 = dx**2 + dy**2 + dz**2
964:                 if(r2 .le. Rcut2)then954:                 if(r2 .le. Rcut2)then
965:                  rm2 = 1/r2955:                  rm2 = 1/r2
966:                  rm14 = rm2**7956:                  rm14 = rm2**7
967:                  energy=energy+STT*rm2**6+SC957:                  energy=energy+STT*rm2**6+SC
968:                  f_over_r=-STT*12.0*rm14958:                  f_over_r=-STT*12.0*rm14
969:                  RD1=sqrt(r2)-NCswitch959:                  RD1=sqrt(r2)-NCswitch
1287: 1277: 
1288: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1278: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1289: !* SBMPR computes the forces due to position restraints               *1279: !* SBMPR computes the forces due to position restraints               *
1290: !**********************************************************************1280: !**********************************************************************
1291: 1281: 
1292:       subroutine SBMPR(x,y,z,grad, energy, natoms, 1282:       subroutine SBMPR(x,y,z,grad, energy, natoms, 
1293:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)1283:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
1294: 1284: 
1295:       USE KEY1285:       USE KEY
1296:       implicit NONE1286:       implicit NONE
1297:       integer I, J, NATOMS,I1,I21287:       integer I, J, NATOMS
1298:       INTEGER SBMPRN, SBMPRI(NATOMS)1288:       INTEGER SBMPRN, SBMPRI(NATOMS)
1299:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)1289:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
1300:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),1290:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
1301:      Q grad(3*NATOMS), energy,ET1291:      Q grad(3*NATOMS), energy,ET
1302:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K61292:       DOUBLE PRECISION DX,DY,DZ,K1,K2,K3,K4,K5,K6
1303: 1293: 
1304: 1294: 
1305:       if(SBMPRN .eq. 0) RETURN1295:       if(SBMPRN .eq. 0) RETURN
1306: 1296: 
1307: !$OMP PARALLEL PRIVATE(I,J,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:energy,grad)1297: !$OMP PARALLEL PRIVATE(I,J,DX,DY,DZ,K1,K2,K3,K4,K5,K6)REDUCTION(+:energy,grad)
1308: !$OMP DO1298: !$OMP DO
1309:               do I=1,SBMPRN1299:               do I=1,SBMPRN
1310:            J=SBMPRI(I)1300:            J=SBMPRI(I)
1311:            I1=(I-1)*61301:            DX=X(J)-SBMPRX(I,1)
1312:            I2=(I-1)*31302:            DY=Y(J)-SBMPRX(I,2)
1313:            DX=X(J)-SBMPRX(I2+1)1303:            DZ=Z(J)-SBMPRX(I,3)
1314:            DY=Y(J)-SBMPRX(I2+2)1304:            K1=SBMPRK(I,1)
1315:            DZ=Z(J)-SBMPRX(I2+3)1305:            K2=SBMPRK(I,2)
1316:            K1=SBMPRK(I1+1)1306:            K3=SBMPRK(I,3)
1317:            K2=SBMPRK(I1+2)1307:            K4=SBMPRK(I,4)
1318:            K3=SBMPRK(I1+3)1308:            K5=SBMPRK(I,5)
1319:            K4=SBMPRK(I1+4)1309:            K6=SBMPRK(I,6)
1320:            K5=SBMPRK(I1+5) 
1321:            K6=SBMPRK(I1+6) 
1322:            1310:            
1323: 1311: 
1324:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +1312:            ENERGY = ENERGY +0.5*(K1*DX**2+K2*DY**2+K3*DZ**2) +
1325:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ1313:      Q     K4*DX*DY+K5*DX*DZ+K6*DY*DZ
1326: 1314: 
1327:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ 1315:            grad(J*3-2) = grad(J*3-2) + K1*DX + K4*DY + K5*DZ 
1328:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ1316:            grad(J*3-1) = grad(J*3-1) + K2*DY + K4*DX + K6*DZ
1329:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY1317:            grad(J*3)   = grad(J*3)   + K3*DZ + K5*DX + K6*DY
1330:         enddo1318:         enddo
1331: !$OMP END DO1319: !$OMP END DO


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0