hdiff output

r30690/SBM.f 2016-07-06 15:38:25.381520581 +0100 r30689/SBM.f 2016-07-06 15:38:25.713525066 +0100
  1:       SUBROUTINE SBM(qo,NATOMS,grad,energy,GTEST,STEST)  1:       SUBROUTINE SBM(qo,NATOMS,grad,energy,GTEST,STEST)
  2:       USE KEY  2:       USE KEY
  3:       implicit NONE  3:       implicit NONE
  4:       INTEGER NATOMS,i,j  4:       INTEGER NATOMS,i,j
  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)  5:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS)
  6:       DOUBLE PRECISION ENERGY  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=12000) 11:         parameter(NSBMMAX=30000)
 12:       INTEGER MAXSEP,MAXSEPSYS 12:       INTEGER MAXSEP,MAXSEPSYS
 13:       PARAMETER (MAXSEP=40) 13:       PARAMETER (MAXSEP=40)
 14:       INTEGER MAXATOMTYPES 14:       INTEGER MAXATOMTYPES
 15:       INTEGER ATOMTYPES 15:       INTEGER ATOMTYPES
 16:       PARAMETER (MAXATOMTYPES=10) 16:       PARAMETER (MAXATOMTYPES=10)
 17:       DOUBLE PRECISION STT,SA, SB, SC 
 18:       DIMENSION STT(MAXATOMTYPES*MAXATOMTYPES),ATOMTYPES(NSBMMAX) 17:       DIMENSION STT(MAXATOMTYPES*MAXATOMTYPES),ATOMTYPES(NSBMMAX)
  18:       DOUBLE PRECISION SA, SB, SC
 19:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES) 19:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
 20:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES) 20:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
 21:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES) 21:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
 22:       INTEGER NNCINC 22:       LOGICAL NNCINC
 23:       DIMENSION NNCINC(NSBMMAX*MAXSEP) 23:       DIMENSION NNCINC(NSBMMAX*MAXSEP)
 24:       DOUBLE PRECISION  Rb(5*NSBMMAX),bK(5*NSBMMAX),ANTC(2*NSBMMAX), 24:       DOUBLE PRECISION  Rb(2*NSBMMAX), bK(2*NSBMMAX), ANTC(2*NSBMMAX),
 25:      Q Tk(2*NSBMMAX),PK(3*NSBMMAX),PHISBM(3*NSBMMAX),CONCOEF(NSBMMAX*4*6), 25:      Q Tk(2*NSBMMAX), PK(3*NSBMMAX), PHISBM(3*NSBMMAX), CONCOEF(NSBMMAX*5*6),
 26:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 26:      Q NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
 27:       INTEGER  Ib1(NSBMMAX*5), Ib2(NSBMMAX*5), IT(NSBMMAX*2), JT(NSBMMAX*2), 27:       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), 28:      Q KT(NSBMMAX*2),IP(NSBMMAX*3), JP(NSBMMAX*3), KP(NSBMMAX*3),
 29:      Q LP(NSBMMAX*3), IC(NSBMMAX*4), JC(NSBMMAX*4), 29:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5),
 30:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX) 30:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX)
 31:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NEXCLUDE,NEXCLUDEELEC 31:       INTEGER MAXEXCLUSIONS,MAXEXCLUSIONSELEC,NEXCLUDE,NEXCLUDEELEC
 32:       PARAMETER (MAXEXCLUSIONS=20) 32:       PARAMETER (MAXEXCLUSIONS=20)
 33:       PARAMETER (MAXEXCLUSIONSELEC=5) 33:       PARAMETER (MAXEXCLUSIONSELEC=5)
 34:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS), 34:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS),
 35:      Q        NNEXL2(NSBMMAX*MAXEXCLUSIONS), 35:      Q NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC),NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC)
 36:      Q        NNEXL1ELEC(NSBMMAX*MAXEXCLUSIONSELEC), 
 37:      Q        NNEXL2ELEC(NSBMMAX*MAXEXCLUSIONSELEC) 
 38:  36: 
 39:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3) 37:       DOUBLE PRECISION SBMPRK(NSBMMAX*6),SBMPRX(NSBMMAX*3)
 40:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX) 38:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX)
 41:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 39:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 42:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 40:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
 43:      Q SBMPRK,SBMPRX,CONCOEF,SA,SB,SC 41:      Q SBMPRK,SBMPRX,CONCOEF,SA,SB,SC
 44:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC, 42:         common /int/ PHITYPE,Ib1,Ib2,IT,JT,KT,IP,JP,KP,LP,IC,JC,
 45:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1, 43:      Q NBA,NTA,NPA,NC,CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,
 46:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,MAXSEPSYS, 44:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,MAXSEPSYS,
 47:      Q ATOMTYPES 45:      Q ATOMTYPES
  46:         common /logical/ NNCINC
 48:  47: 
 49:         if(NATOMS.gt. NSBMMAX)then 48:         if(NATOMS.gt. NSBMMAX)then
 50:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 49:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 51:         STOP 50:         STOP
 52:         endif 51:         endif
 53:  52: 
 54:        if(.NOT.CALLED)then 53:        if(.NOT.CALLED)then
 55:         NEXCLUDE=MAXEXCLUSIONS 54:         NEXCLUDE=MAXEXCLUSIONS
 56:         NEXCLUDEELEC=MAXEXCLUSIONS 55:         NEXCLUDEELEC=MAXEXCLUSIONS
 57:         MAXSEPSYS=0 56:         MAXSEPSYS=0
 58:       DO I=1,NATOMS 
 59:         DO J=1,MAXSEP 
 60:           NNCINC((I-1)*MAXSEP+J)=0 
 61:         enddo 
 62:       enddo 
 63:  
 64:         write(*,*) 57:         write(*,*)
 65:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 58:         write(*,*)  'Calculations will use a Structure-based SMOG model:'
 66:         write(*,*)  '   For more information on SMOG models, see:' 59:         write(*,*)  '   For more information on SMOG models, see:'
 67:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 60:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.'
 68:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 61:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.'
 69:         write(*,*) 62:         write(*,*)
 70:          63:         
 71:  64: 
 72:        call SBMinit(NATOMS,ATOMTYPES,MAXATOMTYPES,Ib1,Ib2,Rb,bK,IT,JT, 65:        call SBMinit(NATOMS,ATOMTYPES,MAXATOMTYPES,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,
 73:      Q KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1, 66:      Q PK,PHITYPE,PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE,
 74:      Q NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP, 67:      Q NNEXL1ELEC,NNEXL2ELEC,NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,
 75:      Q MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,SBMCHARGEON,NCswitch,NCcut, 68:      Q SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,SA,SB,SC,CONTACTTYPE,
 76:      Q STT,SA,SB,SC,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN, 69:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 77:      Q SBMPRI,SBMPRK,SBMPRX) 70:         if(NBA .gt. NSBMMAX*2 .or. NTA .gt. NSBMMAX*2 .or.
 78:         if(NBA .gt. NSBMMAX*5 .or. NTA .gt. NSBMMAX*2 .or. 71:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 5*NSBMMAX)then
 79:      Q  NPA .gt. NSBMMAX*3 .or. NC .gt. 4*NSBMMAX)then 
 80:         write(*,*) 'ERROR: arrays not large large enough for the system' 72:         write(*,*) 'ERROR: arrays not large large enough for the system'
 81:         write(*,*) '   entries=',NBA, NTA,NPA,NC  73:         write(*,*) '   entries=',NBA, NTA,NPA,NC 
 82:         stop 74:         stop
 83:         endif  75:         endif 
 84:         CALLED=.TRUE. 76:         CALLED=.TRUE.
 85:         endIF 77:         endIF
 86:  78: 
 87:  79: 
 88: ! call the energy routine 80: ! call the energy routine
 89:       call calc_energy_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,energy, 81:       call calc_energy_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,energy,
 90:      Q Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC, 82:      Q Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,PHISBM,IC,
 91:      Q JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC, 83:      Q JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUDE,NNEXL1ELEC,NNEXL2ELEC,
 92:      Q NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE, 84:      Q NEXCLUDEELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,SBMCHARGE,
 93:      Q SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,KAPPA, 85:      Q SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,
 94:      Q DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 86:      Q DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 95:  
 96:  
 97:       IF (STEST) THEN 87:       IF (STEST) THEN
 98:          PRINT '(A)','ERROR - second derivatives not available' 88:          PRINT '(A)','ERROR - second derivatives not available'
 99:          STOP 89:          STOP
100:       ENDIF 90:       ENDIF
101:       return 91:       return
102:       end 92:       end
103:  93: 
104:  94: 
105: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 95: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
106: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 96: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
122:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 112:      Q KT1,IT2,JT2,KT2,IP1,JP1,KP1,LP1,IP2,JP2,KP2,LP2,nBA1,nTA1,nPA1, 
123:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,113:      Q nBA2,nTA2,nPA2,ind1,ind2,ANt,MDT1, MDT2, cl1,cl2,tempi,
124:      Q MAXEXCLUSIONS,MAXEXCLUSIONSNEW,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,114:      Q MAXEXCLUSIONS,MAXEXCLUSIONSNEW,MAXEXCLUSIONSELEC,NT,NTE,MAXSEP,
125:      Q MAXSEPSYS,ATOMTYPES(NATOMS),ATYPE,POS,MAXATOMTYPES,NATOMTYPES115:      Q MAXSEPSYS,ATOMTYPES(NATOMS),ATYPE,POS,MAXATOMTYPES,NATOMTYPES
126:       DOUBLE PRECISION :: Rb(2*NATOMS), bK(2*NATOMS), ANTC(2*NATOMS), 116:       DOUBLE PRECISION :: Rb(2*NATOMS), bK(2*NATOMS), ANTC(2*NATOMS), 
127:      Q Tk(2*NATOMS), PK(3*NATOMS), PHISBM(3*NATOMS),117:      Q Tk(2*NATOMS), PK(3*NATOMS), PHISBM(3*NATOMS),
128:      Q Sigma, EpsC, NCswitch,NCcut,TMPSTT(MAXATOMTYPES),118:      Q Sigma, EpsC, NCswitch,NCcut,TMPSTT(MAXATOMTYPES),
129:      Q STT(MAXATOMTYPES**2),SBMCHARGE(NATOMS)119:      Q STT(MAXATOMTYPES**2),SBMCHARGE(NATOMS)
130:       DOUBLE PRECISION RDIFF,ALPHA120:       DOUBLE PRECISION RDIFF,ALPHA
131:       DOUBLE PRECISION CONCOEF(NC*6)121:       DOUBLE PRECISION CONCOEF(NC*6)
132:       INTEGER  Ib1(5*NATOMS), Ib2(5*NATOMS), IT(2*NATOMS), JT(2*NATOMS), 122:       INTEGER  Ib1(2*NATOMS), Ib2(2*NATOMS), IT(2*NATOMS), JT(2*NATOMS), 
133:      Q KT(2*NATOMS),IP(3*NATOMS),JP(3*NATOMS),KP(3*NATOMS),123:      Q KT(2*NATOMS),IP(3*NATOMS),JP(3*NATOMS),KP(3*NATOMS),
134:      Q LP(3*NATOMS),IC(NATOMS*4),JC(NATOMS*4), 124:      Q LP(3*NATOMS),IC(NATOMS*5),JC(NATOMS*5), 
135:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE125:      Q NBA,NTA,NPA,NC,SBMCHARGEON(NATOMS),BONDTYPE
136:       INTEGER PHITYPE(3*NATOMS),M126:       INTEGER PHITYPE(3*NATOMS),M
137:       integer AA,BB,ANTEMP127:       integer AA,BB,ANTEMP
138:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)128:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)
139:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)129:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)
140:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)130:       INTEGER NNEXL1ELEC(NATOMS*MAXEXCLUSIONSELEC)
141:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)131:       INTEGER NNEXL2ELEC(NATOMS*MAXEXCLUSIONSELEC)
142:       DOUBLE PRECISION dx,dy,dz132:       DOUBLE PRECISION dx,dy,dz
143:       double precision PI,ARG1,ARG2,ARG3,ARG4133:       double precision PI,ARG1,ARG2,ARG3,ARG4
144:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut134:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
145:       integer CONTACTTYPE(NATOMS*5),TMPARG135:       integer CONTACTTYPE(NATOMS*5),TMPARG
146:       DOUBLE PRECISION SA, SB, SC136:       DOUBLE PRECISION SA, SB, SC
147:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)137:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
148:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)138:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
149:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)139:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
150: 140: 
151:       character TMPCHAR141:       character TMPCHAR
152:       integer TMPINT,NUMCHARGES,nexc,I1,I2142:       integer TMPINT,NUMCHARGES,nexc,I1,I2
153:       double precision TMPREAL,concentration143:       double precision TMPREAL,concentration
154:       LOGICAL TARR(NATOMS)144:       LOGICAL NNCINC(NATOMS*MAXSEP),TARR(NATOMS)
155:       INTEGER, DIMENSION(NATOMS*MAXSEP) ::  NNCINC145: 
156:       INTEGER SBMPRN, SBMPRI(NATOMS)146:       INTEGER SBMPRN, SBMPRI(NATOMS)
157:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)147:       DOUBLE PRECISION SBMPRK(NATOMS*6),SBMPRX(NATOMS*3)
158:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS148:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS
159:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)149:       DIMENSION EXCLUSIONS(NATOMS*MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)
160: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 150: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
161: !$OMP PARALLEL151: !$OMP PARALLEL
162: !$    NTHREADS = OMP_GET_NUM_THREADS()152: !$    NTHREADS = OMP_GET_NUM_THREADS()
163: !$      TID = OMP_GET_THREAD_NUM()153: !$      TID = OMP_GET_THREAD_NUM()
164: !$    if(TID .eq. 0)then 154: !$    if(TID .eq. 0)then 
165: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS155: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS
166: !$    endif156: !$    endif
167: !$OMP END PARALLEL157: !$OMP END PARALLEL
168: 158: 
 159:       DO I=1,NATOMS
 160:         DO J=1,MAXSEP
 161:           NNCINC((I-1)*MAXSEP+J)=.TRUE.
 162:         enddo
 163:       enddo
169: 164: 
170:       pi = 3.14159265358979323846264338327950288419716939937510165:       pi = 3.14159265358979323846264338327950288419716939937510
171:       MaxCon=NATOMS*5166:       MaxCon=NATOMS*5
172:       do i=1,NATOMS167:       do i=1,NATOMS
173:         NUMOFEXCLUSIONS(I)=0168:         NUMOFEXCLUSIONS(I)=0
174:       enddo169:       enddo
175: 170: 
176: ! These lines read in the parameters.171: ! These lines read in the parameters.
177:         open(30, file='SBM.INP', status='old', access='sequential')172:         open(30, file='SBM.INP', status='old', access='sequential')
178:         write(*,*) 173:         write(*,*) 
413: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in408: ! INCLUDEEXCLUSION helps keep track of which exclusions to include in
414: ! the model409: ! the model
415: !**************************************************************************410: !**************************************************************************
416: 411: 
417:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,412:       SUBROUTINE INCLUDEEXCLUSIONS(NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,
418:      Q NNCINC,MAXSEP,MAXSEPSYS)413:      Q NNCINC,MAXSEP,MAXSEPSYS)
419:       IMPLICIT NONE414:       IMPLICIT NONE
420:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,415:       INTEGER NATOMS,ATOM1,ATOM2,EXCLIST,NEXCLUDE,MAXEXCLUDE,MAXSEP,
421:      Q MAXSEPSYS416:      Q MAXSEPSYS
422:       DIMENSION EXCLIST(NATOMS*MAXEXCLUDE),NEXCLUDE(NATOMS)417:       DIMENSION EXCLIST(NATOMS*MAXEXCLUDE),NEXCLUDE(NATOMS)
423:       INTEGER, DIMENSION(NATOMS*MAXSEP) ::  NNCINC418:       LOGICAL  NNCINC
424:       INTEGER I,N,M,AA,BB,AB,POS419:       DIMENSION NNCINC(NATOMS*MAXSEP)
 420:       INTEGER I,N,M,AA,BB,AB
425:       LOGICAL INCLUDED421:       LOGICAL INCLUDED
 422: 
426: ! If the atoms are within MAXSEP of each other, then keep track via423: ! If the atoms are within MAXSEP of each other, then keep track via
427: ! NNCINC424: ! NNCINC
428:       AA=min(ATOM1,ATOM2)425:       AA=min(ATOM1,ATOM2)
429:       BB=max(ATOM1,ATOM2)426:       BB=max(ATOM1,ATOM2)
430:       AB=BB-AA427:       AB=BB-AA
431:       POS=(AA-1)*MAXSEP+AB 
432:       IF(AB .le. MAXSEP)THEN428:       IF(AB .le. MAXSEP)THEN
433:         NNCINC(POS)=1429:         NNCINC((AA-1)*MAXSEP+AB)=.FALSE.
434:         IF(AB .gt. MAXSEPSYS)THEN430:         IF(AB .gt. MAXSEPSYS)THEN
435:           MAXSEPSYS=AB431:           MAXSEPSYS=AB
436:         ENDIF432:         ENDIF
437:       ELSE433:       ELSE
438:         INCLUDED = .FALSE.434:         INCLUDED = .FALSE.
439:         N=MAX(ATOM1,ATOM2)435:         N=MAX(ATOM1,ATOM2)
440:         M=MIN(ATOM1,ATOM2)436:         M=MIN(ATOM1,ATOM2)
441:         DO I =1,NEXCLUDE(M)437:         DO I =1,NEXCLUDE(M)
442:           IF(EXCLIST((M-1)*MAXEXCLUDE+I) .eq. N)THEN438:           IF(EXCLIST((M-1)*MAXEXCLUDE+I) .eq. N)THEN
443:             INCLUDED = .TRUE.439:             INCLUDED = .TRUE.
445:         ENDDO 441:         ENDDO 
446:   442:   
447:         if(.NOT. INCLUDED)THEN443:         if(.NOT. INCLUDED)THEN
448:           if(NEXCLUDE(M) .EQ. MAXEXCLUDE)THEN444:           if(NEXCLUDE(M) .EQ. MAXEXCLUDE)THEN
449:             write(*,*) 'ERROR: TOO MANY EXCLUSIONS WITH ATOM',M445:             write(*,*) 'ERROR: TOO MANY EXCLUSIONS WITH ATOM',M
450:           ENDIF446:           ENDIF
451:           NEXCLUDE(M)=NEXCLUDE(M)+1447:           NEXCLUDE(M)=NEXCLUDE(M)+1
452:           EXCLIST((M-1)*MAXEXCLUDE+NEXCLUDE(M))=N448:           EXCLIST((M-1)*MAXEXCLUDE+NEXCLUDE(M))=N
453:         ENDIF449:         ENDIF
454:       ENDIF450:       ENDIF
455:  
456:       END451:       END
457: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^452: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
458: 453: 
459: 454: 
460: C455: C
461: C Calculate the Forces and energies456: C Calculate the Forces and energies
462: C457: C
463:       subroutine CALC_ENERGY_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,458:       subroutine CALC_ENERGY_SBM(qo,natoms,ATOMTYPES,MAXATOMTYPES,GRAD,
464:      Q energy,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,459:      Q energy,Ib1,Ib2,Rb,bK,IT,JT,KT,ANTC,Tk,IP,JP,KP,LP,PK,PHITYPE,
465:      Q PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,460:      Q PHISBM,IC,JC,CONCOEF,NNEXL1,NNEXL2,NEXCLUSIONS,NNNEXL1ELEC,
466:      Q NNEXL2ELEC,NEXCLUSIONSELEC,NNCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,461:      Q NNEXL2ELEC,NEXCLUSIONSELEC,NCINC,MAXSEP,MAXSEPSYS,NBA,NTA,NPA,NC,
467:      Q SBMCHARGE,SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,462:      Q SBMCHARGE,SBMCHARGEON,NCswitch,STT,SA,SB,SC,NCcut,CONTACTTYPE,PREFACTOR,
468:      Q KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)463:      Q KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
469: 464: 
470:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,465:       INTEGER I,J,NATOMS,NBA,NTA,NPA,NC,NEXCLUSIONS,MAXSEP,
471:      Q NEXCLUSIONSELEC,ATOMTYPES(NATOMS),MAXATOMTYPES466:      Q NEXCLUSIONSELEC,ATOMTYPES(NATOMS),MAXATOMTYPES
472: 467: 
473:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY468:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
474:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)469:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
475:       INTEGER NNCINC470:       LOGICAL NNCINC
476:       DIMENSION NNCINC(NATOMS*MAXSEP)471:       DIMENSION NNCINC(NATOMS*MAXSEP)
477:       DOUBLE PRECISION SA, SB, SC472:       DOUBLE PRECISION SA, SB, SC
478:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)473:       DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
479:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)474:       DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
480:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)475:       DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
481: 476: 
482:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 477:         DOUBLE PRECISION  Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
483:      Q PHISBM(NPA),NCswitch,NCcut,STT(MAXATOMTYPES**2), SBMCHARGE(NATOMS)478:      Q PHISBM(NPA),NCswitch,NCcut,STT(MAXATOMTYPES**2), SBMCHARGE(NATOMS)
484:         DOUBLE PRECISION  CONCOEF(NC*6)479:         DOUBLE PRECISION  CONCOEF(NC*6)
485:         INTEGER Ib1(NBA),Ib2(NBA),IT(NTA),JT(NTA),KT(NTA),IP(NPA), 480:         INTEGER Ib1(NBA),Ib2(NBA),IT(NTA),JT(NTA),KT(NTA),IP(NPA), 
507:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,502:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
508:      Q PHITYPE,PHISBM,NPA)503:      Q PHITYPE,PHISBM,NPA)
509:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC,CONCOEF,504:       call SBMContacts(x,y,z,grad, energy,natoms,IC,JC,CONCOEF,
510:      Q NC,CONTACTTYPE)505:      Q NC,CONTACTTYPE)
511:       call SBMNonContacts(x,y,z,grad,energy,natoms,ATOMTYPES,506:       call SBMNonContacts(x,y,z,grad,energy,natoms,ATOMTYPES,
512:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,507:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,
513:      Q NCswitch,NCcut,STT,SA,SB,SC)508:      Q NCswitch,NCcut,STT,SA,SB,SC)
514:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,509:       call SBMDHELEC(x,y,z,grad,energy,natoms,NNEXL1ELEC,NNEXL2ELEC,
515:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)510:      Q NEXCLUSIONSELEC,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
516:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)511:       call SBMPR(x,y,z,grad, energy, natoms,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 512: 
517:       end513:       end
518: 514: 
519: 515: 
520: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<516: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
521: !* SBMBonds  computes the hookean force and energy between chosen atoms *517: !* SBMBonds  computes the hookean force and energy between chosen atoms *
522: !***********************************************************************518: !***********************************************************************
523: 519: 
524:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)520:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)
525:       USE KEY521:       USE KEY
526:       implicit NONE522:       implicit NONE
936: !* SBMNonContacts computes the forces due to non native contacts      *932: !* SBMNonContacts computes the forces due to non native contacts      *
937: !**********************************************************************933: !**********************************************************************
938: 934: 
939:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,ATOMTYPES,935:         subroutine SBMnoncontacts(x,y,z,grad,energy,NATOMS,ATOMTYPES,
940:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,936:      Q MAXATOMTYPES,NNEXL1,NNEXL2,NEXCLUSIONS,NNCINC,MAXSEP,MAXSEPSYS,
941:      Q NCswitch,NCcut,STT,SA,SB,SC)937:      Q NCswitch,NCcut,STT,SA,SB,SC)
942:         USE KEY938:         USE KEY
943:         implicit NONE939:         implicit NONE
944:         integer I,N,J,AN,NATOMS,ATOMTYPES(NATOMS),MAXATOMTYPES940:         integer I,N,J,AN,NATOMS,ATOMTYPES(NATOMS),MAXATOMTYPES
945: 941: 
 942: 
946:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),943:         DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
947:      Q  grad(3*NATOMS), energy944:      Q  grad(3*NATOMS), energy,STT(MAXATOMTYPES**2),ET
948:        DOUBLE PRECISION STT,SA,SB,SC,SAT,SBT,SCT,STTT945:        DOUBLE PRECISION SA, SB, SC,SAT,SBT,SCT,STTT
949:        DIMENSION STT(MAXATOMTYPES*MAXATOMTYPES) 
950:        DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)946:        DIMENSION SA(MAXATOMTYPES*MAXATOMTYPES)
951:        DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)947:        DIMENSION SB(MAXATOMTYPES*MAXATOMTYPES)
952:        DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)948:        DIMENSION SC(MAXATOMTYPES*MAXATOMTYPES)
953: 949: 
954: 950: 
955:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj,POS951:         integer C1,C2,DINDEX,ii,jj,kk,k,l,iii,jjj,POS
956:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 952:         DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
957:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,MAXSEPSYS953:         INTEGER NTHREADS,NEXCLUSIONS,MAXSEP,MAXSEPSYS
958:         INTEGER NNCINC954:         LOGICAL NNCINC
959:         DIMENSION NNCINC(NATOMS*MAXSEP)955:         DIMENSION NNCINC(NATOMS*MAXSEP)
960: !$    INTEGER  OMP_GET_NUM_THREADS,956: !$    INTEGER  OMP_GET_NUM_THREADS,
961: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS957: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
962: 958: 
963:         INTEGER NNEXL1(NEXCLUSIONS)959:         INTEGER NNEXL1(NEXCLUSIONS)
964:         INTEGER NNEXL2(NEXCLUSIONS)960:         INTEGER NNEXL2(NEXCLUSIONS)
965:         DOUBLE PRECISION dx,dy,dz961:         DOUBLE PRECISION dx,dy,dz
966:         integer tempN, alpha962:         integer tempN, alpha
967:         double precision Rdiff,Vfunc,Ffunc963:         double precision Rdiff,Vfunc,Ffunc
968:         double precision Rcut2,Rswitch2964:         double precision Rcut2,Rswitch2
978:         integer MaxGridX,MaxGridY,MaxGridZ974:         integer MaxGridX,MaxGridY,MaxGridZ
979:         double precision gridsize,RD1975:         double precision gridsize,RD1
980:         double precision minX,minY,minZ,maxX,maxY,maxZ976:         double precision minX,minY,minZ,maxX,maxY,maxZ
981:         integer Xgrid,Ygrid,Zgrid,TID977:         integer Xgrid,Ygrid,Zgrid,TID
982:         Rdiff=NCcut-NCswitch978:         Rdiff=NCcut-NCswitch
983:         alpha=12979:         alpha=12
984:         GRIDSIZE=NCcut*1.01980:         GRIDSIZE=NCcut*1.01
985:         Rcut2=NCcut**2981:         Rcut2=NCcut**2
986:         Rswitch2=NCswitch**2982:         Rswitch2=NCswitch**2
987: 983: 
 984: !        SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
 985: !     Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
 986: 
 987: !        SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)
 988: 
 989: !        SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)
988: ! grid the system990: ! grid the system
989:         minX=10000000991:         minX=10000000
990:         minY=10000000992:         minY=10000000
991:         minZ=10000000993:         minZ=10000000
992:         994:         
993:         maxX=-10000000995:         maxX=-10000000
994:         maxY=-10000000996:         maxY=-10000000
995:         maxZ=-10000000997:         maxZ=-10000000
996: 998: 
997:         do i=1,NATOMS999:         do i=1,NATOMS
1070:           Zgrid=int((Z(I)-minZ)/gridsize)+11072:           Zgrid=int((Z(I)-minZ)/gridsize)+1
1071:           do ii=XGRID-1,XGRID+11073:           do ii=XGRID-1,XGRID+1
1072:            do jj=YGRID-1,YGRID+11074:            do jj=YGRID-1,YGRID+1
1073:             do kk=ZGRID-1,ZGRID+11075:             do kk=ZGRID-1,ZGRID+1
1074:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.1076:              if(ii .ge. 1 .and. jj .ge. 1 .and. kk .ge. 1 .and.
1075:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then1077:      Q            ii .le. MAXGRIDX .and. jj .le. MAXGRIDY .and. kk .le.MAXGRIDZ)then
1076:               do jjj=1,Ngrid(ii,jj,kk)1078:               do jjj=1,Ngrid(ii,jj,kk)
1077:                C2=grid(ii,jj,kk,jjj)1079:                C2=grid(ii,jj,kk,jjj)
1078:                DINDEX=C2-I1080:                DINDEX=C2-I
1079:                if(DINDEX .GT. 0)THEN1081:                if(DINDEX .GT. 0)THEN
1080:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX) .eq. 0) )then1082:                 if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS .AND. NNCINC((I-1)*MAXSEP+DINDEX)  ) )then
1081:                 !if(DINDEX .GT. MAXSEPSYS .OR. (DINDEX .LE. MAXSEPSYS )  )then 
1082:                  dx = X(I) - X(C2)1083:                  dx = X(I) - X(C2)
1083:                  dy = Y(I) - Y(C2)1084:                  dy = Y(I) - Y(C2)
1084:                  dz = Z(I) - Z(C2)1085:                  dz = Z(I) - Z(C2)
1085:                  r2 = dx**2 + dy**2 + dz**21086:                  r2 = dx**2 + dy**2 + dz**2
1086:                  if(r2 .le. Rcut2)then1087:                  if(r2 .le. Rcut2)then
1087:                   POS=(ATOMTYPES(I)-1)*MAXATOMTYPES+ATOMTYPES(C2)1088:                   POS=(ATOMTYPES(I)-1)*MAXATOMTYPES+ATOMTYPES(C2)
1088:                   STTT=STT(POS)1089:                   STTT=STT(POS)
1089:                   SCT=SC(POS)1090:                   SCT=SC(POS)
1090:                   rm2 = 1/r21091:                   rm2 = 1/r2
1091:                   rm14 = rm2**71092:                   rm14 = rm2**7


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0