hdiff output

r30655/SBM.f 2016-07-06 15:35:42.403317173 +0100 r30654/SBM.f 2016-07-06 15:35:42.747321825 +0100
 24:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP, 24:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP,
 25:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC,NNCinc, 25:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC,NNCinc,
 26:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI 26:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI
 27:         if(NATOMS.gt. NSBMMAX)then 27:         if(NATOMS.gt. NSBMMAX)then
 28:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX' 28:         write(*,*) 'TOO MANY ATOMS FOR SBM, change NSBMMAX'
 29:         STOP 29:         STOP
 30:         endif 30:         endif
 31:  31: 
 32:        if(.NOT.CALLED)then 32:        if(.NOT.CALLED)then
 33:         write(*,*) 33:         write(*,*)
 34:         write(*,*)  'Calculations will use a Structure-based SMOG model:' 34:         write(*,*)  'Calculations will use a Structure-based SMOG model, described in:'
 35:         write(*,*)  '   Software: Noel, et al. PLoS Comput Biol 12, e1004794, 2016.' 35:         write(*,*)  'Whitford, et al. Prot. Struct. Func. Bioinfo. 75, 430-441, 2009.'
 36:         write(*,*)  '   Model: Whitford, et al. Prot Struct Func Bioinfo 75, 430-441, 2009.' 
 37:         write(*,*) 36:         write(*,*)
 38:          37:         
 39:  38: 
 40:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP,  39:        call SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, IP, 
 41:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, Sigma,  40:      Q JP, KP, LP, PK, PHITYPE,PHISBM,IC, JC, Sigma, 
 42:      Q EpsC, NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT, 41:      Q EpsC, NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,NCswitch,NCcut,STT,
 43:      Q NSBMMAX,CONTACTTYPE, 42:      Q NSBMMAX,CONTACTTYPE,
 44:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX) 43:      Q PREFACTOR,KAPPA,DHswitch,DHcut,SBMPRN,SBMPRI,SBMPRK,SBMPRX)
 45:         if(NBA .gt. NSBMMAX .or. NTA .gt. NSBMMAX .or. 44:         if(NBA .gt. NSBMMAX .or. NTA .gt. NSBMMAX .or.
 46:      Q  NPA .gt. NSBMMAX .or. NC .gt. 5*NSBMMAX)then 45:      Q  NPA .gt. NSBMMAX .or. NC .gt. 5*NSBMMAX)then
 51:         CALLED=.TRUE. 50:         CALLED=.TRUE.
 52:         endIF 51:         endIF
 53: ! call the energy routine 52: ! call the energy routine
 54:  53: 
 55:       call calc_energy_SBM(qo,natoms, GRAD, energy, Ib1, Ib2, 54:       call calc_energy_SBM(qo,natoms, GRAD, energy, Ib1, Ib2,
 56:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK, 55:      Q Rb, bK, IT, JT, KT, ANTC, Tk, IP, JP, KP, LP, PK,
 57:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC,  56:      Q PHITYPE,PHISBM,IC, JC, Sigma, EpsC, 
 58:      Q  NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON, 57:      Q  NNCinc,NBA, NTA, NPA, NC, SBMCHARGE,SBMCHARGEON,
 59:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut, 58:      Q NCswitch,STT,NCcut,CONTACTTYPE,PREFACTOR,KAPPA,DHswitch,DHcut,
 60:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX) 59:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
  60: 
 61:       IF (STEST) THEN 61:       IF (STEST) THEN
 62:          PRINT '(A)','ERROR - second derivatives not available' 62:          PRINT '(A)','ERROR - second derivatives not available'
 63:          STOP 63:          STOP
 64:       ENDIF 64:       ENDIF
 65:       return 65:       return
 66:       end 66:       end
 67:  67: 
 68:  68: 
  69: 
  70: 
  71: 
 69: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 72: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 70: !* SBMinit() reads the atom positions from file.  If 1 is selected for * 73: !* SBMinit() reads the atom positions from file.  If 1 is selected for *
 71: !* startt then the velocities are assigned, otherwise, they are read   * 74: !* startt then the velocities are assigned, otherwise, they are read   *
 72: !* by selecting 2, or generated by selecting 3                         * 75: !* by selecting 2, or generated by selecting 3                         *
 73: !*********************************************************************** 76: !***********************************************************************
 74:  77: 
 75:       subroutine SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk, 78:       subroutine SBMinit(NATOMS,Ib1, Ib2,Rb, bK, IT, JT, KT, ANTC, Tk,
 76:      Q  IP, JP, KP, LP, PK,PHITYPE,PHISBM,IC, JC, Sigma,  79:      Q  IP, JP, KP, LP, PK,PHITYPE,PHISBM,IC, JC, Sigma, 
 77:      Q EpsC, NNCinc, NBA,  80:      Q EpsC, NNCinc, NBA, 
 78:      Q NTA, NPA, NC, SBMCHARGE, SBMCHARGEON,NCswitch,NCcut,STT,NSBMMAX, 81:      Q NTA, NPA, NC, SBMCHARGE, SBMCHARGEON,NCswitch,NCcut,STT,NSBMMAX,
106:         logical TEMPARRAY109:         logical TEMPARRAY
107:         dimension TEMPARRAY(NATOMS,NATOMS)110:         dimension TEMPARRAY(NATOMS,NATOMS)
108:         integer CONTACTTYPE111:         integer CONTACTTYPE
109:         character TMPCHAR112:         character TMPCHAR
110:         integer TMPINT,NUMCHARGES,nexc,I1,I2113:         integer TMPINT,NUMCHARGES,nexc,I1,I2
111:         double precision TMPREAL,concentration114:         double precision TMPREAL,concentration
112: 115: 
113:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)116:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)
114:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)117:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
115: 118: 
116: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 119: 
117: !$OMP PARALLEL120: 
118: !$    NTHREADS = OMP_GET_NUM_THREADS() 
119: !$      TID = OMP_GET_THREAD_NUM() 
120: !$    if(TID .eq. 0)then  
121: !$      write(*,*) 'OMP enabled. Number of threads:', NTHREADS 
122: !$    endif 
123: !$OMP END PARALLEL 
124:       pi = 3.14159265358979323846264338327950288419716939937510121:       pi = 3.14159265358979323846264338327950288419716939937510
125: 122: 
126:       NNCmax = NATOMS*NATOMS123:         NNCmax = NATOMS*NATOMS
127:       MaxCon=NATOMS*5124:         MaxCon=NATOMS*5
128: 125: 
129:       do i=1,NATOMS126:         do i=1,NATOMS
130:         do j=1,NATOMS127:                 do j=1,NATOMS
131:           TEMPARRAY(i,j)=.TRUE.128:                 TEMPARRAY(i,j)=.TRUE.
 129:                 enddo
132:         enddo130:         enddo
133:       enddo 
134: 131: 
135: 132: 
136: 133: 
137: ! These lines read in the parameters.134: ! These lines read in the parameters.
138:         open(30, file='SBM.INP', status='old', access='sequential')135:         open(30, file='SBM.INP', status='old', access='sequential')
139:         write(*,*) 136:         write(*,*) 
140:         write(*,*) 'Reading the forcefield from SBM.INP'137:         write(*,*) 'Reading the forcefield from SBM.INP'
141:         write(*,*) 138:         write(*,*) 
142:         read(30,*)139:         read(30,*)
143:         read(30,*)140:         read(30,*)
144:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut141:         read(30,*) PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut
145:         read(30,*)142:         read(30,*)
146:         read(30,*) RSig, Reps, NCswitch, NCcut143:         read(30,*) RSig, Reps,NCswitch,NCcut
147:         write(*,*) 'RSig, Reps,NCswitch,NCcut'144:         write(*,*) 'RSig, Reps,NCswitch,NCcut'
148:         write(*,'(4F10.5)') RSig, Reps, NCswitch, NCcut145:         write(*,*) RSig, Reps,NCswitch,NCcut
149:         write(*,*) 'Reading electrostatic parameters'146:         write(*,*) 'Reading electrostatic parameters'
150:         write(*,'(5F10.5)') PREFACTOR,DC,CONCENTRATION,DHswitch,DHcut147:         write(*,*) PREFACTOR,DC,KAPPA,DHswitch,DHcut
151:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is is units148:         ! CONCENTRATION IS THE MONOVALENT ION CONCENTRATION kappa is is units
152:         ! A^-1149:         ! A^-1
153:         KAPPA=0.32*sqrt(CONCENTRATION)150:         KAPPA=0.32*sqrt(CONCENTRATION)
154:         PREFACTOR=PREFACTOR/DC151:         PREFACTOR=PREFACTOR/DC
155:         read(30,*) ANtemp152:         read(30,*) ANtemp
156:         write(*,*) ANtemp, 'atoms'153:         write(*,*) ANtemp, ' atoms'
157:         NUMCHARGES=1154:         NUMCHARGES=1
158:         do i=1, ANtemp155:         do i=1, ANtemp
159:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)156:           read(30,*) TMPINT,TMPCHAR,TMPINT,TMPCHAR,TMPCHAR,SBMCHARGE(i),ATMASS(I)
160:           if(SBMCHARGE(i) .ne. 0)then157:           if(SBMCHARGE(i) .ne. 0)then
161:              NUMCHARGES=NUMCHARGES+1158:              NUMCHARGES=NUMCHARGES+1
162:              SBMCHARGEON(NUMCHARGES)=i159:              SBMCHARGEON(NUMCHARGES)=i
163:           endif 160:           endif 
164:         end do161:         end do
165:         SBMCHARGEON(1)=NUMCHARGES-1162:         SBMCHARGEON(1)=NUMCHARGES-1
166:         163:         
253: 250: 
254:                 STT=reps*rsig**12251:                 STT=reps*rsig**12
255: 252: 
256:         do I=1,NATOMS253:         do I=1,NATOMS
257:           TARR(I)=0254:           TARR(I)=0
258:         enddo255:         enddo
259: 256: 
260: ! read in position restraints257: ! read in position restraints
261:        read(30,*) 258:        read(30,*) 
262:        read(30,*) SBMPRN259:        read(30,*) SBMPRN
263:        write(*,*) SBMPRN, 'position restraints'260:        write(*,*) SBMPRN, ' position restraints'
264:        do I=1,SBMPRN261:        do I=1,SBMPRN
265:             read(30,*) SBMPRI(I),SBMPRK(I,1),SBMPRK(I,2),SBMPRK(I,3),262:             read(30,*) SBMPRI(I),SBMPRK(I,1),SBMPRK(I,2),SBMPRK(I,3),
266:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6),263:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6),
267:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)264:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)
268:             if(TARR(SBMPRI(I)) .ne. 0)then265:             if(TARR(SBMPRI(I)) .ne. 0)then
269:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)266:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
270: !                call abort267: !                call abort
271:             endif268:             endif
272:             TARR(SBMPRI(I))=1269:             TARR(SBMPRI(I))=1
273:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  270:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  
311:         y(i) = qo(j+2)308:         y(i) = qo(j+2)
312:         z(i) = qo(j+3)309:         z(i) = qo(j+3)
313:         grad(j+1) = 0.0310:         grad(j+1) = 0.0
314:         grad(j+2) = 0.0311:         grad(j+2) = 0.0
315:         grad(j+3) = 0.0312:         grad(j+3) = 0.0
316:       enddo313:       enddo
317: 314: 
318:       energy = 0.0315:       energy = 0.0
319:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)316:       call SBMbonds(x,y,z,grad, energy, natoms,Ib1, Ib2,Rb, bK,NBA)
320:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)317:       call SBMangl(x,y,z,grad, energy, natoms,IT,JT,KT,ANTC,Tk,NTA)
321:       call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,318:         call SBMDihedral(x,y,z,grad, energy, natoms,IP,JP,KP,LP,PK,
322:      Q PHITYPE,PHISBM,NPA)319:      Q PHITYPE,PHISBM,NPA)
323:       call SBMContacts(x,y,z,grad, energy, natoms, IC, JC, 320:         call SBMContacts(x,y,z,grad, energy, natoms, IC, JC, 
324:      Q Sigma, EpsC, NC,CONTACTTYPE)321:      Q Sigma, EpsC, NC,CONTACTTYPE)
325:       call SBMNonContacts(x,y,z,grad, energy, natoms, 322:         call SBMNonContacts(x,y,z,grad, energy, natoms, 
326:      Q NNCinc,NCswitch,NCcut,STT)323:      Q NNCinc,NCswitch,NCcut,STT)
327:       call SBMDHELEC(x,y,z,grad, energy, natoms, 324:         call SBMDHELEC(x,y,z,grad, energy, natoms, 
328:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)325:      Q NNCinc,PREFACTOR,KAPPA,DHswitch,DHcut,SBMCHARGE,SBMCHARGEON)
329:       call SBMPR(x,y,z,grad, energy, natoms, 326:         call SBMPR(x,y,z,grad, energy, natoms, 
330:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)327:      Q SBMPRN,SBMPRI,SBMPRK,SBMPRX)
331: 328: 
332:       end329:       end
333: 330: 
334: 331: 
335: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<332: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
336: !* SBMBonds  computes the hookean force and energy between chosen atoms *333: !* SBMBonds  computes the hookean force and energy between chosen atoms *
337: !***********************************************************************334: !***********************************************************************
338: 335: 
339:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)336:       subroutine SBMBonds(x,y,z,grad,energy, natoms,Ib1, Ib2,Rb, bK,NBA)
508: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<505: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
509: !* SBMdihedral computes the dihedral angles and the forces due to them *506: !* SBMdihedral computes the dihedral angles and the forces due to them *
510: !**********************************************************************507: !**********************************************************************
511: 508: 
512:       SUBROUTINE SBMdihedral(x,y,z,grad, energy, NATOMS,IP,JP,KP,LP,PK,509:       SUBROUTINE SBMdihedral(x,y,z,grad, energy, NATOMS,IP,JP,KP,LP,PK,
513:      Q PHITYPE,PHISBM,NPA)510:      Q PHITYPE,PHISBM,NPA)
514:       USE KEY511:       USE KEY
515:       implicit NONE512:       implicit NONE
516:       integer I, N, J, NATOMS, NPA, II513:       integer I, N, J, NATOMS, NPA, II
517:       DOUBLE PRECISION x(NATOMS),y(NATOMS),z(NATOMS),514:       DOUBLE PRECISION x(NATOMS),y(NATOMS),z(NATOMS),
518:      Q grad(3*NATOMS),gradOMP(3*NATOMS),energy,ETEMP515:      Q grad(3*NATOMS),energy
519: 516: 
520: 517: 
521:       DOUBLE PRECISION PK(NPA),PHISBM(NPA)518:       DOUBLE PRECISION PK(NPA),PHISBM(NPA)
522:       INTEGER IP(NPA), JP(NPA), KP(NPA), LP(NPA) 519:       INTEGER IP(NPA), JP(NPA), KP(NPA), LP(NPA) 
523:       INTEGER PHITYPE(NPA)520:       INTEGER PHITYPE(NPA)
524: 521: 
525:       double precision lfac522:       double precision lfac
526:       integer I3, J3, K3, L3,C1523:       integer I3, J3, K3, L3
527:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,524:       double precision  XIJ,YIJ,ZIJ,XKJ,YKJ,
528:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,525:      + ZKJ,XKL,YKL,ZKL,RIJ, RKJ,RKL,DX,DY,
529:      + DZ, GX,GY,GZ,CT,CPHI,526:      + DZ, GX,GY,GZ,CT,CPHI,
530:      + SPHI,Z1, Z2,FXI,FYI,FZI,527:      + SPHI,Z1, Z2,FXI,FYI,FZI,
531:      + FXJ,FYJ,FZJ, FXK,FYK,FZK,528:      + FXJ,FYJ,FZJ, FXK,FYK,FZK,
532:      + FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,ftem,CT0,CT1,AP0,AP1,529:      + FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,ftem,CT0,CT1,AP0,AP1,
533:      + Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,530:      + Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ,
534:      + S,HGoverG,FGoverG,A1,A3531:      + S,HGoverG,FGoverG,A1,A3
535: 532: 
536: 533:       DIMENSION XIJ(NPA),YIJ(NPA),ZIJ(NPA),RIJ(NPA),XKJ(NPA),
 534:      + YKJ(NPA),
 535:      + ZKJ(NPA),RKJ(NPA),XKL(NPA),YKL(NPA),ZKL(NPA),RKL(NPA),DX(NPA),
 536:      + DY(NPA),
 537:      + DZ(NPA), GX(NPA),GY(NPA),GZ(NPA),CT(NPA),CPHI(NPA),
 538:      + SPHI(NPA),Z1(NPA), Z2(NPA),FXI(NPA),FYI(NPA),FZI(NPA),
 539:      + FXJ(NPA),FYJ(NPA),FZJ(NPA), FXK(NPA),FYK(NPA),FZK(NPA),
 540:      + FXL(NPA),FYL(NPA),FZL(NPA),DF(NPA),FGoverG(NPA),HGoverG(NPA)
537: C541: C
538:       double precision  TM24,TM06,tenm3,zero,one,two,four,six,twelve542:       double precision  TM24,TM06,tenm3,zero,one,two,four,six,twelve
539: 543: 
540:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,544:       DOUBLE PRECISION TT1, TT2, TT3, TT4, TT1X,TT1Y,TT1Z,TT2X,TT2Y,
541:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z545:      + TT2Z, TT3X, TT3Y, TT3Z, TT4X, TT4Y, TT4Z
542: 546: 
543:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/547:       DATA TM24,TM06,tenm3/1.0d-24,1.0d-06,1.0d-03/
544:       data zero,one,two,four,six,twelve/0.d0,1.d0,2.d0,4.d0,6.d0,12.d0/548:       data zero,one,two,four,six,twelve/0.d0,1.d0,2.d0,4.d0,6.d0,12.d0/
545: 549: 
546:       double precision pi550:       double precision pi
547:       pi = 3.14159265358979323846264338327950288419716939937510551:       pi = 3.14159265358979323846264338327950288419716939937510
548:       do i=1,NATOMS 
549:          gradOMP(i*3)=0 
550:          gradOMP(i*3-1)=0 
551:          gradOMP(i*3-2)=0 
552:       enddo 
553: 552: 
554:       etemp=0 
555: !$OMP PARALLEL PRIVATE(I,I3,J3,K3,L3,II,XIJ,YIJ,ZIJ,XKJ,YKJ,  
556: !$OMP& ZKJ,XKL,YKL,ZKL,RIJ,RKJ,RKL,DX,DY,  
557: !$OMP& DZ, GX,GY,GZ,CT,CPHI, 
558: !$OMP& SPHI,Z1, Z2,FXI,FYI,FZI, 
559: !$OMP& FXJ,FYJ,FZJ, FXK,FYK,FZK, 
560: !$OMP& FXL,FYL,FZL,DF,Z10,Z20,Z12,Z11,Z22,CT0,CT1,AP0,AP1, 
561: !$OMP& Dums,DFLIM, DF1, DF0, DR1, DR2,DR3,DR4,DR5,DR6,DRX,DRY,DRZ, 
562: !$OMP& S,HGoverG,FGoverG,A1,A3,TT1,TT2,TT3,TT4,TT1X,TT1Y,TT1Z,TT2X, 
563: !$OMP& TT2Y,TT2Z,TT3X,TT3Y,TT3Z,TT4X,TT4Y,TT4Z) REDUCTION(+:etemp,gradOMP) 
564: !! 
565: !$OMP DO 
566:           DO II = 1,nPA553:           DO II = 1,nPA
567: 554: 
568:             I3 = IP(II)555:             I3 = IP(II)
569:             J3 = JP(II)556:             J3 = JP(II)
570:             K3 = KP(II)557:             K3 = KP(II)
571:             L3 = LP(II)558:             L3 = LP(II)
572:             XIJ = X(I3)-X(J3)559: 
573:             YIJ = Y(I3)-Y(J3)560:  
574:             ZIJ = Z(I3)-Z(J3)561: 
575:             XKJ = X(K3)-X(J3)562:             XIJ(II) = X(I3)-X(J3)
576:             YKJ = Y(K3)-Y(J3)563:             YIJ(II) = Y(I3)-Y(J3)
577:             ZKJ = Z(K3)-Z(J3)564:             ZIJ(II) = Z(I3)-Z(J3)
578:             RKJ = sqrt(XKJ**2+YKJ**2+ZKJ**2)565:             XKJ(II) = X(K3)-X(J3)
579:             XKL = X(K3)-X(L3)566:             YKJ(II) = Y(K3)-Y(J3)
580:             YKL = Y(K3)-Y(L3)567:             ZKJ(II) = Z(K3)-Z(J3)
581:             ZKL = Z(K3)-Z(L3)                                  568:             RKJ(II) = sqrt(XKJ(II)**2+YKJ(II)**2+ZKJ(II)**2)
582:             FGoverG=-(XIJ*XKJ+YIJ*YKJ+ZIJ*ZKJ)/RKJ569:             XKL(II) = X(K3)-X(L3)
583:             HGoverG=(XKL*XKJ+YKL*YKJ+ ZKL*ZKJ)/RKJ570:             YKL(II) = Y(K3)-Y(L3)
 571:             ZKL(II) = Z(K3)-Z(L3)                                  
 572: 
 573: 
 574:             FGoverG(II)=-(XIJ(II)*XKJ(II)+YIJ(II)*YKJ(II)+
 575:      Q   ZIJ(II)*ZKJ(II))/RKJ(II)
 576:             HGoverG(II)=(XKL(II)*XKJ(II)+YKL(II)*YKJ(II)+
 577:      Q   ZKL(II)*ZKJ(II))/RKJ(II)
 578:           END DO
584: C DX is the M vector and G is the N vector579: C DX is the M vector and G is the N vector
585:             DX = YIJ*ZKJ-ZIJ*YKJ580:           DO II = 1,nPA
586:             DY = ZIJ*XKJ-XIJ*ZKJ581:             DX(II) = YIJ(II)*ZKJ(II)-ZIJ(II)*YKJ(II)
587:             DZ = XIJ*YKJ-YIJ*XKJ582:             DY(II) = ZIJ(II)*XKJ(II)-XIJ(II)*ZKJ(II)
588:             GX = ZKJ*YKL-YKJ*ZKL583:             DZ(II) = XIJ(II)*YKJ(II)-YIJ(II)*XKJ(II)
589:             GY = XKJ*ZKL-ZKJ*XKL584:             GX(II) = ZKJ(II)*YKL(II)-YKJ(II)*ZKL(II)
590:             GZ = YKJ*XKL-XKJ*YKL585:             GY(II) = XKJ(II)*ZKL(II)-ZKJ(II)*XKL(II)
591:             FXI = SQRT(DX*DX+DY*DY+DZ*DZ)586:             GZ(II) = YKJ(II)*XKL(II)-XKJ(II)*YKL(II)
592:             FYI = SQRT(GX*GX+GY*GY+GZ*GZ)587:           END DO
593:             CT = DX*GX+DY*GY+DZ*GZ588: C
594:             z10 = 1.0/FXI589: 
595:             z20 = 1.0/FYI590: 
 591: ! so far so good
 592: 
 593: 
 594:           DO II = 1,nPA
 595:             FXI(II) = SQRT(DX(II)*DX(II)
 596:      Q                    +DY(II)*DY(II)
 597:      Q                    +DZ(II)*DZ(II))
 598:             FYI(II) = SQRT(GX(II)*GX(II)
 599:      Q                    +GY(II)*GY(II)
 600:      Q                    +GZ(II)*GZ(II))
 601:             CT(II) = DX(II)*GX(II)+DY(II)*GY(II)+DZ(II)*GZ(II)
 602:           END DO
 603: 
 604:          DO II = 1,nPA
 605:             z10 = 1.0/FXI(II)
 606:             z20 = 1.0/FYI(II)
596:             Z12 = Z10*Z20607:             Z12 = Z10*Z20
597:             Z1 = Z10608:             Z1(II) = Z10
598:             Z2 = Z20609:             Z2(II) = Z20
599:             CT0 = MIN(1.0,CT*Z12)610:             ftem = zero
600:             CT1 = MAX(-1.0,CT0)611:             CT0 = MIN(one,CT(II)*Z12)
601:             S = XKJ*(DZ*GY-DY*GZ)+YKJ*(DX*GZ-DZ*GX)+ZKJ*(DY*GX-DX*GY)612:             CT1 = MAX(-one,CT0)
 613:             S = XKJ(II)*(DZ(II)*GY(II)-DY(II)*GZ(II))+
 614:      Q          YKJ(II)*(DX(II)*GZ(II)-DZ(II)*GX(II))+
 615:      Q          ZKJ(II)*(DY(II)*GX(II)-DX(II)*GY(II))
602:             AP0 = ACOS(CT1)616:             AP0 = ACOS(CT1)
603:             AP1 = PI-SIGN(AP0,S)617:             AP1 = PI-SIGN(AP0,S)
604: 618: 
605: 619:             CT(II) = AP1
606:             CT = AP1620:             CPHI(II) = COS(AP1)
607:             CPHI = COS(AP1)621:             SPHI(II) = SIN(AP1)
608:             SPHI = SIN(AP1)622:          END DO
609: 623: 
610: 624: 
 625:         DO II = 1,nPA
611: !            CT0 = CT(II)626: !            CT0 = CT(II)
612: ! Here is the energy part627: ! Here is the energy part
613:           A1=CT-PHISBM(II)628:           A1=CT(II)-PHISBM(II)
614:           A3=A1*3629:           A3=A1*3
615: 630: 
616:         if(PHITYPE(II) .eq. 1)then631:         if(PHITYPE(II) .eq. 1)then
617:           etemp =  etemp + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))632:           Energy =  Energy + PK(II)*(3.0/2.0-cos(A1)-0.5*cos(A3))
618: ! dE/dPHI633: ! dE/dPHI
619:             DF=PK(II)*(sin(A1)+1.5*sin(A3))634:             DF(II)=PK(II)*(sin(A1)+1.5*sin(A3))
620:         elseif(PHITYPE(II) .eq. 2)then635:         elseif(PHITYPE(II) .eq. 2)then
621:           if(A1 .gt. PI)then636:           if(A1 .gt. PI)then
622:                 A1=A1-2*PI637:                 A1=A1-2*PI
623:           elseif(A1 .lt. -PI)then638:           elseif(A1 .lt. -PI)then
624:                 A1=A1+2*PI639:                 A1=A1+2*PI
625:           endif640:           endif
626: 641: 
627:           etemp =  etemp + 0.5*PK(II)*A1**2642:           Energy =  Energy + 0.5*PK(II)*A1**2
628: ! dE/dPHI643: ! dE/dPHI
629:             DF=PK(II)*A1644:             DF(II)=PK(II)*A1
630: 645: 
631:         else646:         else
632: 647: 
633:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)648:         write(*,*) 'unrecognized dihedral type', PHITYPE(II)
634:         STOP649:         STOP
635:         endif650:         endif
 651:         END DO
 652: 
 653: 
636: 654: 
637: ! insert the new 655: ! insert the new 
638: 656: 
639: ! now, do dPhi/dX657:        ! now, do dPhi/dX
 658: 
 659:          DO II = 1,nPA
640: 660: 
641: ! |G|/|A|**2 661: ! |G|/|A|**2 
642:             TT1 = Z1*Z1*RKJ*DF662:             TT1 = Z1(II)*Z1(II)*RKJ(II)*DF(II)
643: ! FG/(A**2*|G|)663: ! FG/(A**2*|G|)
644:             TT2 = FGoverG*Z1*Z1*DF664:             TT2 = FGoverG(II)*Z1(II)*Z1(II)*DF(II)
645: ! HG/(B**2*|G|)665: ! HG/(B**2*|G|)
646:             TT3 = HGoverG*Z2*Z2*DF666:             TT3 = HGoverG(II)*Z2(II)*Z2(II)*DF(II)
647: ! |G|/|B|**2 667: ! |G|/|B|**2 
648:             TT4 = Z2*Z2*RKJ*DF668:             TT4 = Z2(II)*Z2(II)*RKJ(II)*DF(II)
 669: 
 670: 
649: ! note: negatives are flipped from paper because A=-DX671: ! note: negatives are flipped from paper because A=-DX
650:         TT1X=TT1*DX672:         TT1X=TT1*DX(II)
651:         TT1Y=TT1*DY673:         TT1Y=TT1*DY(II)
652:         TT1Z=TT1*DZ674:         TT1Z=TT1*DZ(II)
653: 675: 
654:         TT2X=TT2*DX676:         TT2X=TT2*DX(II)
655:         TT2Y=TT2*DY677:         TT2Y=TT2*DY(II)
656:         TT2Z=TT2*DZ678:         TT2Z=TT2*DZ(II)
657: 679: 
658: 680: 
659:         TT3X=TT3*GX681:         TT3X=TT3*GX(II)
660:         TT3Y=TT3*GY682:         TT3Y=TT3*GY(II)
661:         TT3Z=TT3*GZ683:         TT3Z=TT3*GZ(II)
662:  
663:  
664:         TT4X=TT4*GX 
665:         TT4Y=TT4*GY 
666:         TT4Z=TT4*GZ 
667:  
668:             gradOMP(I3*3-2) =  gradOMP(I3*3-2)  + TT1X   
669:             gradOMP(I3*3-1) =  gradOMP(I3*3-1)  + TT1Y 
670:             gradOMP(I3*3)   =  gradOMP(I3*3)    + TT1Z 
671:             gradOMP(J3*3-2) =  gradOMP(J3*3-2)  - TT1X - TT2X - TT3X 
672:             gradOMP(J3*3-1) =  gradOMP(J3*3-1)  - TT1Y - TT2Y - TT3Y 
673:             gradOMP(J3*3)   =  gradOMP(J3*3)    - TT1Z - TT2Z - TT3Z 
674:             gradOMP(K3*3-2) =  gradOMP(K3*3-2)  + TT2X + TT3X - TT4X 
675:             gradOMP(K3*3-1) =  gradOMP(K3*3-1)  + TT2Y + TT3Y - TT4Y 
676:             gradOMP(K3*3)   =  gradOMP(K3*3)    + TT2Z + TT3Z - TT4Z 
677:             gradOMP(L3*3-2) =  gradOMP(L3*3-2)  + TT4X 
678:             gradOMP(L3*3-1) =  gradOMP(L3*3-1)  + TT4Y 
679:             gradOMP(L3*3)   =  gradOMP(L3*3)    + TT4Z 
680:         END DO 
681: !$OMP END DO  
682: !$OMP END PARALLEL 
683: 684: 
684:         ENERGY=ENERGY+etemp 
685:         do C1=1,NATOMS 
686:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2)  
687:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1)  
688:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)    
689:         enddo 
690: 685: 
 686:         TT4X=TT4*GX(II)
 687:         TT4Y=TT4*GY(II)
 688:         TT4Z=TT4*GZ(II)
 689: 
 690:             I3 = IP(II)
 691:             J3 = JP(II)
 692:             K3 = KP(II)
 693:             L3 = LP(II)
 694: 
 695: 
 696:             grad(I3*3-2) =  grad(I3*3-2)  + TT1X  
 697:             grad(I3*3-1) =  grad(I3*3-1)  + TT1Y
 698:             grad(I3*3)   =  grad(I3*3)    + TT1Z
 699:             grad(J3*3-2) =  grad(J3*3-2)  - TT1X - TT2X - TT3X
 700:             grad(J3*3-1) =  grad(J3*3-1)  - TT1Y - TT2Y - TT3Y
 701:             grad(J3*3)   =  grad(J3*3)    - TT1Z - TT2Z - TT3Z
 702:             grad(K3*3-2) =  grad(K3*3-2)  + TT2X + TT3X - TT4X
 703:             grad(K3*3-1) =  grad(K3*3-1)  + TT2Y + TT3Y - TT4Y
 704:             grad(K3*3)   =  grad(K3*3)    + TT2Z + TT3Z - TT4Z
 705:             grad(L3*3-2) =  grad(L3*3-2)  + TT4X
 706:             grad(L3*3-1) =  grad(L3*3-1)  + TT4Y
 707:             grad(L3*3)   =  grad(L3*3)    + TT4Z
 708: 
 709:           END DO
691: 710: 
692:           END711:           END
693: 712: 
694: 713: 
695: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^714: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^END of SBMDihedral^^^^^^^^^^^^^^^^^^^^^^^^^^^
696: 715: 
697: 716: 
698: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<717: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
699: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *718: !* SBMCONTACTS: computes the force on all atoms due to contacts via a   *
700: !* 10-12 or 6-12 potential                                              *719: !* 10-12 or 6-12 potential                                              *
801:       implicit NONE820:       implicit NONE
802:       integer I, N, J, AN, NATOMS821:       integer I, N, J, AN, NATOMS
803: 822: 
804: 823: 
805:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),824:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS),
806:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC825:      Q grad(3*NATOMS), gradOMP(3*NATOMS), energy,STT,SA,SB,SC
807: 826: 
808:       integer C1, C2, ii,jj,kk,k,l,iii,jjj827:       integer C1, C2, ii,jj,kk,k,l,iii,jjj
809:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 828:       DOUBLE PRECISION  r2, rm2, rm14, f_over_r, NCswitch,NCcut 
810: 829: 
811: !$    INTEGER NTHREADS, OMP_GET_NUM_THREADS,830:         INTEGER NTHREADS,  OMP_GET_NUM_THREADS,
812: !$   Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS831:      Q  OMP_GET_THREAD_NUM,OMP_NUM_THREADS
813: 832: 
814:         INTEGER NNCinc(NATOMS,NATOMS)833:         INTEGER NNCinc(NATOMS,NATOMS)
815:         DOUBLE PRECISION dx,dy,dz834:         DOUBLE PRECISION dx,dy,dz
816:         integer tempN, alpha835:         integer tempN, alpha
817:         double precision Rdiff,Vfunc,Ffunc836:         double precision Rdiff,Vfunc,Ffunc
818:         double precision Rcut2,Rswitch2837:         double precision Rcut2,Rswitch2
819:         ! Ngrid is the number of atoms in that grid point838:         ! Ngrid is the number of atoms in that grid point
820:         ! grid is the array of atoms in each grid839:         ! grid is the array of atoms in each grid
821:         integer Ngrid,grid,maxgrid,maxpergrid840:         integer Ngrid,grid,maxgrid,maxpergrid
822:         ! number of atoms per grid, max841:         ! number of atoms per grid, max
838:         Rcut2=NCcut**2857:         Rcut2=NCcut**2
839:         Rswitch2=NCswitch**2858:         Rswitch2=NCswitch**2
840: 859: 
841:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 860:         SB=-1.0/Rdiff**3*( 2*alpha*STT/NCcut**(alpha+1)  + 
842:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )861:      Q    (alpha)*(alpha+1)*STT*Rdiff/NCcut**(alpha+2)   )
843: 862: 
844:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)863:         SA=-(alpha*(alpha+1)*STT/NCcut**(alpha+2)+3*SB*Rdiff**2)/(2*Rdiff)
845: 864: 
846:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)865:         SC=-(STT/NCcut**alpha+SA/3.0*Rdiff**3+SB/4.0*Rdiff**4)
847: 866: 
 867: 
 868: 
848: ! grid the system869: ! grid the system
849:         minX=10000000870:         minX=10000000
850:         minY=10000000871:         minY=10000000
851:         minZ=10000000872:         minZ=10000000
852:         873:         
853:         maxX=-10000000874:         maxX=-10000000
854:         maxY=-10000000875:         maxY=-10000000
855:         maxZ=-10000000876:         maxZ=-10000000
856: 877: 
857:         do i=1,NATOMS878:         do i=1,NATOMS
 879:         gradOMP(i*3)=0
 880:         gradOMP(i*3-1)=0
 881:         gradOMP(i*3-2)=0
858: 882: 
859:            gradOMP(i*3)=0 
860:            gradOMP(i*3-1)=0 
861:            gradOMP(i*3-2)=0 
862:            if(X(i) .gt. maxX)then883:            if(X(i) .gt. maxX)then
863:                 maxX=X(i)884:                 maxX=X(i)
864:            elseif(Y(i) .gt. maxY)then885:            elseif(Y(i) .gt. maxY)then
865:                 maxY=Y(i)886:                 maxY=Y(i)
866:            endif887:            endif
867: 888: 
868:            if(Z(i) .gt. maxZ)then889:            if(Z(i) .gt. maxZ)then
869:                 maxZ=Z(i)890:                 maxZ=Z(i)
870:            elseif(X(i) .lt. minX)then891:            elseif(X(i) .lt. minX)then
871:                 minX=X(i)892:                 minX=X(i)
968:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx989:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx
969:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy990:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy
970:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz991:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz
971: 992: 
972:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx993:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx
973:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy994:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy
974:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz995:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz
975: 996: 
976:             endif997:             endif
977: 998: 
 999: !!! END
978: 1000: 
979:            endif1001:            endif
980:            enddo1002:            enddo
981:  1003:  
982:           endif1004:           endif
983: 1005: 
984:             enddo1006:             enddo
985:            enddo1007:            enddo
986:           enddo1008:           enddo
987: 1009: 
988:           enddo1010:           enddo
989: !$OMP END PARALLEL1011: !$OMP END PARALLEL
990:           energy = energy  + etemp1012:                 energy = energy  + etemp
991:         do C1=1,NATOMS1013:         do C1=1,NATOMS
992:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 1014:               grad(C1*3-2) = grad(C1*3-2)+ gradOMP(C1*3-2) 
993:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 1015:               grad(C1*3-1) = grad(C1*3-1)+ gradOMP(C1*3-1) 
994:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   1016:               grad(C1*3)   = grad(C1*3)  + gradOMP(C1*3)   
995:         enddo1017:         enddo
996: 1018: 
997:       END1019:       END
998: 1020: 
999: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^1021: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
1000: 1022: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0