hdiff output

r30678/SBM.f 2016-07-06 15:38:20.425453607 +0100 r30677/SBM.f 2016-07-06 15:38:20.753458029 +0100
 17:      Q CONCOEF2(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut 17:      Q CONCOEF2(NSBMMAX*5),NCswitch,NCcut,SBMCHARGE(NSBMMAX),PREFACTOR,KAPPA,DHswitch,DHcut
 18:       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),
 19:      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),
 20:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5), 20:      Q LP(NSBMMAX*3), IC(NSBMMAX*5), JC(NSBMMAX*5),
 21:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX) 21:      Q PHITYPE(3*NSBMMAX),NBA,NTA,NPA,NC,SBMPRN, SBMPRI(NSBMMAX)
 22:       INTEGER MAXEXCLUSIONS,NEXCLUDE 22:       INTEGER MAXEXCLUSIONS,NEXCLUDE
 23:       PARAMETER (MAXEXCLUSIONS=20) 23:       PARAMETER (MAXEXCLUSIONS=20)
 24:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS) 24:       INTEGER NNEXL1(NSBMMAX*MAXEXCLUSIONS),NNEXL2(NSBMMAX*MAXEXCLUSIONS)
 25:  25: 
 26:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3) 26:       DOUBLE PRECISION SBMPRK(NSBMMAX,6),SBMPRX(NSBMMAX,3)
 27:         integer CONTACTTYPE(NSBMMAX*5),SBMCHARGEON(NSBMMAX) 27:         integer CONTACTTYPE,SBMCHARGEON(NSBMMAX)
 28:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM, 28:         common /double precision/ Rb, bK, ANTC, Tk,PK, PHISBM,
 29:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut, 29:      Q NCswitch,NCcut,STT,SBMCHARGE,PREFACTOR,KAPPA,DHswitch,DHcut,
 30:      Q SBMPRK,SBMPRX,CONCOEF1, CONCOEF2 30:      Q SBMPRK,SBMPRX,CONCOEF1, CONCOEF2
 31:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP, 31:         common /int/ PHITYPE, Ib1, Ib2, IT, JT, KT, IP,
 32:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC, 32:      Q JP, KP, LP,IC, JC, NBA, NTA, NPA, NC,
 33:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,NNEXL2,NEXCLUDE 33:      Q CONTACTTYPE,SBMCHARGEON,SBMPRN,SBMPRI,NNEXL1,NNEXL2,NEXCLUDE
 34:         common /logical/ NNCINC 34:         common /logical/ NNCINC
 35:         MAXSEP=MAXSEP1 35:         MAXSEP=MAXSEP1
 36:  36: 
 37:         if(NATOMS.gt. NSBMMAX)then 37:         if(NATOMS.gt. NSBMMAX)then
106:      Q KT(2*NSBMMAX),IP(3*NSBMMAX), JP(3*NSBMMAX), KP(3*NSBMMAX),106:      Q KT(2*NSBMMAX),IP(3*NSBMMAX), JP(3*NSBMMAX), KP(3*NSBMMAX),
107:      Q LP(3*NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 107:      Q LP(3*NSBMMAX), IC(NSBMMAX*5), JC(NSBMMAX*5), 
108:      Q  NBA, NTA, NPA, NC,SBMCHARGEON(NATOMS)108:      Q  NBA, NTA, NPA, NC,SBMCHARGEON(NATOMS)
109:       INTEGER PHITYPE(3*NSBMMAX)109:       INTEGER PHITYPE(3*NSBMMAX)
110:       integer AA,BB,ANTEMP110:       integer AA,BB,ANTEMP
111:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)111:       INTEGER NNEXL1(NATOMS*MAXEXCLUSIONS)
112:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)112:       INTEGER NNEXL2(NATOMS*MAXEXCLUSIONS)
113:       DOUBLE PRECISION dx,dy,dz113:       DOUBLE PRECISION dx,dy,dz
114:       double precision PI114:       double precision PI
115:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut115:       DOUBLE PRECISION RSig, Reps,PREFACTOR,DC,KAPPA,DHswitch,DHcut
116:       integer CONTACTTYPE(NSBMMAX*5)116:       integer CONTACTTYPE
117:       character TMPCHAR117:       character TMPCHAR
118:       integer TMPINT,NUMCHARGES,nexc,I1,I2118:       integer TMPINT,NUMCHARGES,nexc,I1,I2
119:       double precision TMPREAL,concentration119:       double precision TMPREAL,concentration
120:       LOGICAL NNCINC(NATOMS,MAXSEP)120:       LOGICAL NNCINC(NATOMS,MAXSEP)
121: 121: 
122:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)122:       INTEGER SBMPRN, SBMPRI(NATOMS),TARR(NATOMS)
123:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)123:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
124:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS124:       INTEGER EXCLUSIONS,NUMOFEXCLUSIONS
125:       DIMENSION EXCLUSIONS(NATOMS,MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)125:       DIMENSION EXCLUSIONS(NATOMS,MAXEXCLUSIONS),NUMOFEXCLUSIONS(NATOMS)
126: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 126: !$    INTEGER NTHREADS,OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM,TID 
174:           if(SBMCHARGE(i) .ne. 0)then174:           if(SBMCHARGE(i) .ne. 0)then
175:              NUMCHARGES=NUMCHARGES+1175:              NUMCHARGES=NUMCHARGES+1
176:              SBMCHARGEON(NUMCHARGES)=i176:              SBMCHARGEON(NUMCHARGES)=i
177:           endif 177:           endif 
178:         end do178:         end do
179:         SBMCHARGEON(1)=NUMCHARGES179:         SBMCHARGEON(1)=NUMCHARGES
180:         180:         
181:         read(30,*) 181:         read(30,*) 
182:         read(30,*) NC182:         read(30,*) NC
183:         write(*,*) NC, 'contacts'        183:         write(*,*) NC, 'contacts'        
 184:         read(30,*) CONTACTTYPE
184:           if(NC .gt. MaxCon)then185:           if(NC .gt. MaxCon)then
185:              write(*,*) 'too many contacts'186:              write(*,*) 'too many contacts'
186:              STOP187:              STOP
187:           endif188:           endif
188: 189: 
189:         do i=1, NC190:         do i=1, NC
190:           read(30, *) IC(i), JC(i), CONTACTTYPE(I),Sigma, EpsC191:           read(30, *) IC(i), JC(i), Sigma, EpsC
191:           !! Add more contact types192:           !! Add more contact types
192:           CALL INCLUDEEXCLUSIONS(NATOMS,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)193:           CALL INCLUDEEXCLUSIONS(NATOMS,IC(I),JC(I),EXCLUSIONS,NUMOFEXCLUSIONS,MAXEXCLUSIONS,NNCINC,MAXSEP,MAXSEPTEMP)
193:           if(CONTACTTYPE(I) .eq. 1)then194:           if(CONTACTTYPE .eq. 1)then
194:             ConCoef1(I)=EPSC*SIGMA**12-REPS*RSIG**12195:             ConCoef1(I)=EPSC*SIGMA**12-REPS*RSIG**12
195:             ConCoef2(I)=2*EPSC*SIGMA**6196:             ConCoef2(I)=2*EPSC*SIGMA**6
196:           elseif(CONTACTTYPE(I) .eq. 2)then197:           elseif(CONTACTTYPE .eq. 2)then
197:             ConCoef1(I)=5*EPSC*SIGMA**12-REPS*RSIG**12198:             ConCoef1(I)=5*EPSC*SIGMA**12-REPS*RSIG**12
198:             ConCoef2(I)=6*EPSC*SIGMA**10199:             ConCoef2(I)=6*EPSC*SIGMA**10
199:           else200:           else
200:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.'201:             write(*,*) 'ERROR: Only contacttypes 1 and 2 supported.'
201:           endif202:           endif
202:         end do203:         end do
203: ! make a function to check and add204: ! make a function to check and add
204: 205: 
205:           read(30,*)206:           read(30,*)
206:           read(30,*) nBA207:           read(30,*) nBA
294: ! read in position restraints295: ! read in position restraints
295:        read(30,*) 296:        read(30,*) 
296:        read(30,*) SBMPRN297:        read(30,*) SBMPRN
297:        write(*,*) SBMPRN, 'position restraints'298:        write(*,*) SBMPRN, 'position restraints'
298:        do I=1,SBMPRN299:        do I=1,SBMPRN
299:             read(30,*) SBMPRI(I),SBMPRK(I,1),SBMPRK(I,2),SBMPRK(I,3),300:             read(30,*) SBMPRI(I),SBMPRK(I,1),SBMPRK(I,2),SBMPRK(I,3),
300:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6),301:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6),
301:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)302:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)
302:             if(TARR(SBMPRI(I)) .ne. 0)then303:             if(TARR(SBMPRI(I)) .ne. 0)then
303:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)304:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
304:                 STOP305: !                call abort
305:             endif306:             endif
306:             TARR(SBMPRI(I))=1307:             TARR(SBMPRI(I))=1
307:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  308:         if(SBMPRK(I,4) .ne. 0 .or. SBMPRK(I,5) .ne. 0 .or. SBMPRK(I,6) .ne.0)then  
308:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'309:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'
309:         endif310:         endif
310:        enddo 311:        enddo 
311:         MAXSEP=MAXSEPTEMP312:         MAXSEP=MAXSEPTEMP
312:        close(30)313:        close(30)
313:       END 314:       END 
314: 315: 
375: 376: 
376:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY377:       DOUBLE PRECISION qo(3*NATOMS), grad(3*NATOMS), ENERGY
377:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)378:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS)
378:       LOGICAL NNCINC(NATOMS,MAXSEP)379:       LOGICAL NNCINC(NATOMS,MAXSEP)
379: 380: 
380:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 381:         DOUBLE PRECISION Rb(NBA), bK(NBA), ANTC(NTA), Tk(NTA), PK(NPA), 
381:      Q PHISBM(NPA),CONCOEF1(NC),CONCOEF2(NC),NCswitch,NCcut,STT, 382:      Q PHISBM(NPA),CONCOEF1(NC),CONCOEF2(NC),NCswitch,NCcut,STT, 
382:      Q SBMCHARGE(NATOMS)383:      Q SBMCHARGE(NATOMS)
383:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 384:         INTEGER Ib1(NBA), Ib2(NBA), IT(NTA), JT(NTA), KT(NTA),IP(NPA), 
384:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),SBMCHARGEON(NATOMS),385:      Q JP(NPA), KP(NPA), LP(NPA), IC(NC), JC(NC),SBMCHARGEON(NATOMS),
385:      Q PHITYPE(NPA),CONTACTTYPE(NC)386:      Q PHITYPE(NPA),CONTACTTYPE
386:         INTEGER NNEXL1(NEXCLUSIONS)387:         INTEGER NNEXL1(NEXCLUSIONS)
387:         INTEGER NNEXL2(NEXCLUSIONS)388:         INTEGER NNEXL2(NEXCLUSIONS)
388:       DOUBLE PRECISION dx,dy,dz389:       DOUBLE PRECISION dx,dy,dz
389:       INTEGER SBMPRN, SBMPRI(NATOMS)390:       INTEGER SBMPRN, SBMPRI(NATOMS)
390:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)391:       DOUBLE PRECISION SBMPRK(NATOMS,6),SBMPRX(NATOMS,3)
391:       do i = 1, natoms392:       do i = 1, natoms
392:         j = (i-1)*3393:         j = (i-1)*3
393:         x(i) = qo(j+1)394:         x(i) = qo(j+1)
394:         y(i) = qo(j+2)395:         y(i) = qo(j+2)
395:         z(i) = qo(j+3)396:         z(i) = qo(j+3)
719:       integer I, NATOMS,NC720:       integer I, NATOMS,NC
720: 721: 
721:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 722:       DOUBLE PRECISION x(NATOMS), y(NATOMS), z(NATOMS), 
722:      Q grad(3*NATOMS), energy723:      Q grad(3*NATOMS), energy
723:       DOUBLE PRECISION dx,dy,dz724:       DOUBLE PRECISION dx,dy,dz
724: 725: 
725:       integer C1, C2 726:       integer C1, C2 
726:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 727:       DOUBLE PRECISION  r2, rm2, rm10, rm6,f_over_r,r1 
727: 728: 
728:         DOUBLE PRECISION CONCOEF1(NC), CONCOEF2(NC)729:         DOUBLE PRECISION CONCOEF1(NC), CONCOEF2(NC)
729:         INTEGER IC(NC), JC(NC),CONTACTTYPE(NC)730:         INTEGER IC(NC), JC(NC),CONTACTTYPE
730: 731: 
731: ! type 1 is 6-12732: ! type 1 is 6-12
732: ! type 2 is 12-12733: ! type 2 is 12-12
733: ! implement type 3 as gaussian734: ! implement type 3 as gaussian
734: ! implement type 4 as dual-basin gaussian735: ! implement type 4 as dual-basin gaussian
735: ! question: Should type be part of the interaction, or should they be organized736: ! question: Should type be part of the interaction, or should they be organized
736: ! into separate lists?737: ! into separate lists?
737: 738: 
738: ! type 1 is 6-12 interaction739: ! type 1 is 6-12 interaction
 740:         if(CONTACTTYPE .eq. 1)then
739: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R) REDUCTION(+:ENERGY,grad)741: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM6,F_OVER_R) REDUCTION(+:ENERGY,grad)
740: !$OMP DO742: !$OMP DO
741:           do i=1, NC743:           do i=1, NC
742:             C1 = IC(i)744:             C1 = IC(i)
743:             C2 = JC(i)745:             C2 = JC(i)
744:             dx = X(C1) - X(C2)746:             dx = X(C1) - X(C2)
745:             dy = Y(C1) - Y(C2)747:             dy = Y(C1) - Y(C2)
746:             dz = Z(C1) - Z(C2)748:             dz = Z(C1) - Z(C2)
747:             r2 = dx**2 + dy**2 + dz**2749:             r2 = dx**2 + dy**2 + dz**2
748:             rm2 = 1.0/r2750:             rm2 = 1.0/r2
749: 751:           !  rm2 = rm2*sigma(i)
750:             if(CONTACTTYPE(I) .eq. 1)then752:             rm6 = rm2**3
751:               ! type 1 is 6-12 interaction753:             !ENERGY = ENERGY + epsC(i)*rm6*(rm6-2.0)
752:               rm6 = rm2**3754:             ENERGY = ENERGY + rm6*(CONCOEF1(I)*rm6-CONCOEF2(I))
753:               ENERGY = ENERGY + rm6*(CONCOEF1(I)*rm6-CONCOEF2(I))755:             f_over_r = -12.0*rm6*(CONCOEF1(I)*rm6-CONCOEF2(I)/2.0)*rm2
754:               f_over_r = -12.0*rm6*(CONCOEF1(I)*rm6-CONCOEF2(I)/2.0)*rm2 
755:             elseif(CONTACTTYPE(I) .eq. 2)then 
756:               ! type 2 is 10-12 interaction 
757:               rm10 = rm2**5 
758:               ENERGY = ENERGY + rm10*(CONCOEF1(I)*rm2-CONCOEF2(I)) 
759:               f_over_r = -60.0*rm10*(CONCOEF1(I)*rm2/5.0-CONCOEF2(I)/6.0)*rm2 
760:             else 
761:               WRITE(*,*) 'ERROR: Contact type not recognized:', CONTACTTYPE(I) 
762:               STOP 
763:             endif 
764: 756: 
765:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx757:             grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
766:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy758:             grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
767:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz759:             grad(3*C1)   = grad(3*C1)   + f_over_r * dz
768:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx760:             grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
769:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy761:             grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
770:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz762:             grad(3*C2)   = grad(3*C2)   - f_over_r * dz
771:           enddo763:           enddo
772: !$OMP END DO764: !$OMP END DO
773: !$OMP END PARALLEL765: !$OMP END PARALLEL
 766: ! type 2 is 10-12 interaction
 767:         elseif(CONTACTTYPE .eq. 2)then
 768: !$OMP PARALLEL PRIVATE(I,C1,C2,DX,DY,DZ,R2,RM2,RM10,F_OVER_R) REDUCTION(+:ENERGY,grad)
 769: !$OMP DO
 770:           do i=1, NC
 771:              C1 = IC(i)
 772:              C2 = JC(i)
 773:              dx = X(C1) - X(C2)
 774:              dy = Y(C1) - Y(C2)
 775:              dz = Z(C1) - Z(C2)
 776:              r2 = dx**2 + dy**2 + dz**2
 777:              rm2 = 1.0/r2
 778: !             rm2 = rm2*sigma(i)
 779:              rm10 = rm2**5
 780:              !ENERGY = ENERGY + epsC(i)*rm10*(5*rm2-6.0)
 781:              ENERGY = ENERGY + rm10*(CONCOEF1(I)*rm2-CONCOEF2(I))
 782:              f_over_r = -60.0*rm10*(CONCOEF1(I)*rm2/5.0-CONCOEF2(I)/6.0)*rm2
 783:              grad(3*C1-2) = grad(3*C1-2) + f_over_r * dx
 784:              grad(3*C1-1) = grad(3*C1-1) + f_over_r * dy
 785:              grad(3*C1)   = grad(3*C1)   + f_over_r * dz
 786:              grad(3*C2-2) = grad(3*C2-2) - f_over_r * dx
 787:              grad(3*C2-1) = grad(3*C2-1) - f_over_r * dy
 788:              grad(3*C2)   = grad(3*C2)   - f_over_r * dz
 789:            enddo
 790: !$OMP END DO
 791: !$OMP END PARALLEL
 792: 
 793:         else
 794:         write(*,*) CONTACTTYPE, 'is not a valid contact selection'
 795:         stop
 796:         endif
774: 797: 
775: 798: 
776:       end799:       end
777: 800: 
778: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^801: !^^^^^^^^^^^^^^^^^^^^^^^^^^^^end of SBMContacts^^^^^^^^^^^^^^^^^^^^^^^^^^^
779: 802: 
780: 803: 
781: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<804: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
782: !* SBMNonContacts computes the forces due to non native contacts       *805: !* SBMNonContacts computes the forces due to non native contacts       *
783: !**********************************************************************806: !**********************************************************************
863:         enddo886:         enddo
864: 887: 
865:         maxgridX=int((maxX-minX)/gridsize)+1888:         maxgridX=int((maxX-minX)/gridsize)+1
866:         maxgridY=int((maxY-minY)/gridsize)+1889:         maxgridY=int((maxY-minY)/gridsize)+1
867:         maxgridZ=int((maxZ-minZ)/gridsize)+1890:         maxgridZ=int((maxZ-minZ)/gridsize)+1
868: 891: 
869:         if(maxgridX .ge. maxgrid .or. 892:         if(maxgridX .ge. maxgrid .or. 
870:      Q  maxgridY .ge. maxgrid .or.893:      Q  maxgridY .ge. maxgrid .or.
871:      Q  maxgridZ .ge. maxgrid )then894:      Q  maxgridZ .ge. maxgrid )then
872:         write(*,*) 'ERROR: system got too big for grid searching...'895:         write(*,*) 'ERROR: system got too big for grid searching...'
873:         STOP896: !        call abort
874:         endif897:         endif
875: 898: 
876: 899: 
877: !!  Add a second grid that only includes atoms that are charged900: !!  Add a second grid that only includes atoms that are charged
878: 901: 
879:         do i=1,maxgrid902:         do i=1,maxgrid
880:          do j=1,maxgrid903:          do j=1,maxgrid
881:           do k=1,maxgrid904:           do k=1,maxgrid
882:                 Ngrid(i,j,k)=0905:                 Ngrid(i,j,k)=0
883:           enddo906:           enddo
885:         enddo908:         enddo
886:         do i=1,NATOMS909:         do i=1,NATOMS
887: 910: 
888:                 Xgrid=int((X(i)-minX)/gridsize)+1911:                 Xgrid=int((X(i)-minX)/gridsize)+1
889:                 Ygrid=int((Y(i)-minY)/gridsize)+1912:                 Ygrid=int((Y(i)-minY)/gridsize)+1
890:                 Zgrid=int((Z(i)-minZ)/gridsize)+1913:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
891:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1914:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
892:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then915:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
893:                         write(*,*) 'ERROR: too many atoms in a grid'916:                         write(*,*) 'ERROR: too many atoms in a grid'
894:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid917:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
895:                         STOP918: !                        call abort
896:                 endif919:                 endif
897:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i920:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
898:         enddo921:         enddo
899: 922: 
900: 923: 
901:            tempN=0924:            tempN=0
902: 925: 
903: ! add a second loop that goes over charged atoms926: ! add a second loop that goes over charged atoms
904: 927: 
905: !$OMP PARALLEL928: !$OMP PARALLEL
952:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy975:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
953:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz976:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz
954:                  endif977:                  endif
955:                endif978:                endif
956:               enddo979:               enddo
957:              endif980:              endif
958:             enddo981:             enddo
959:            enddo982:            enddo
960:           enddo983:           enddo
961:          enddo984:          enddo
 985: !$OMP END PARALLEL
962: 986: 
963: !$OMP DO 
964:        DO I=1,NEXCLUSIONS987:        DO I=1,NEXCLUSIONS
965:          C1=NNEXL1(I) 988:          C1=NNEXL1(I) 
966:          C2=NNEXL2(I)989:          C2=NNEXL2(I)
967:          dx = X(C1) - X(C2)990:          dx = X(C1) - X(C2)
968:          dy = Y(C1) - Y(C2)991:          dy = Y(C1) - Y(C2)
969:          dz = Z(C1) - Z(C2)992:          dz = Z(C1) - Z(C2)
970:          r2 = dx**2 + dy**2 + dz**2993:          r2 = dx**2 + dy**2 + dz**2
971:          if(r2 .le. Rcut2)then994:          if(r2 .le. Rcut2)then
972:           rm2 = 1/r2995:           rm2 = 1/r2
973:           rm14 = rm2**7996:           rm14 = rm2**7
984:            write(*,*) 'something went wrong with switching function'1007:            write(*,*) 'something went wrong with switching function'
985:           endif1008:           endif
986:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx1009:           grad(C1*3-2) = grad(C1*3-2) - f_over_r * dx
987:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy1010:           grad(C1*3-1) = grad(C1*3-1) - f_over_r * dy
988:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz1011:           grad(C1*3)   = grad(C1*3)   - f_over_r * dz
989:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx1012:           grad(C2*3-2) = grad(C2*3-2) + f_over_r * dx
990:           grad(C2*3-1) = grad(C2*3-1) + f_over_r * dy1013:           grad(C2*3-1) = grad(C2*3-1) + f_over_r * dy
991:           grad(C2*3)   = grad(C2*3)   + f_over_r * dz1014:           grad(C2*3)   = grad(C2*3)   + f_over_r * dz
992:          endif1015:          endif
993:       ENDDO1016:       ENDDO
994: !$OMP END DO1017: 
995: !$OMP END PARALLEL 
996:       END1018:       END
997: 1019: 
998: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^1020: !^^^^^^^^^^^^^^^^^^^^^^^^^^^End of SBMNonContacts^^^^^^^^^^^^^^^^^^^^^^^^^
999: 1021: 
1000: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<1022: !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1001: ! FUNCTIONS for DH Calculations1023: ! FUNCTIONS for DH Calculations
1002: !*****************************************************1024: !*****************************************************
1003: 1025: 
1004:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)1026:       DOUBLE PRECISION FUNCTION SBMDHV(kappa,r)
1005:       USE KEY1027:       USE KEY
1108:         enddo1130:         enddo
1109: 1131: 
1110:         maxgridX=int((maxX-minX)/gridsize)+11132:         maxgridX=int((maxX-minX)/gridsize)+1
1111:         maxgridY=int((maxY-minY)/gridsize)+11133:         maxgridY=int((maxY-minY)/gridsize)+1
1112:         maxgridZ=int((maxZ-minZ)/gridsize)+11134:         maxgridZ=int((maxZ-minZ)/gridsize)+1
1113: 1135: 
1114:         if(maxgridX .ge. maxgrid .or. 1136:         if(maxgridX .ge. maxgrid .or. 
1115:      Q  maxgridY .ge. maxgrid .or.1137:      Q  maxgridY .ge. maxgrid .or.
1116:      Q  maxgridZ .ge. maxgrid )then1138:      Q  maxgridZ .ge. maxgrid )then
1117:         write(*,*) 'system got too big for grid searching...'1139:         write(*,*) 'system got too big for grid searching...'
1118:         STOP1140: !        call abort
1119:         endif1141:         endif
1120: 1142: 
1121: 1143: 
1122: !!  Add a second grid that only includes atoms that are charged1144: !!  Add a second grid that only includes atoms that are charged
1123: 1145: 
1124:         do i=1,maxgrid1146:         do i=1,maxgrid
1125:          do j=1,maxgrid1147:          do j=1,maxgrid
1126:           do k=1,maxgrid1148:           do k=1,maxgrid
1127:                 Ngrid(i,j,k)=01149:                 Ngrid(i,j,k)=0
1128:           enddo1150:           enddo
1133:         do ii=2,SBMCHARGEON(1)1155:         do ii=2,SBMCHARGEON(1)
1134:                 i=SBMCHARGEON(ii)1156:                 i=SBMCHARGEON(ii)
1135: 1157: 
1136:                 Xgrid=int((X(i)-minX)/gridsize)+11158:                 Xgrid=int((X(i)-minX)/gridsize)+1
1137:                 Ygrid=int((Y(i)-minY)/gridsize)+11159:                 Ygrid=int((Y(i)-minY)/gridsize)+1
1138:                 Zgrid=int((Z(i)-minZ)/gridsize)+11160:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
1139:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+11161:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
1140:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then1162:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
1141:                         write(*,*) 'ERROR: too many atoms in a grid in DH'1163:                         write(*,*) 'ERROR: too many atoms in a grid in DH'
1142:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid1164:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
1143:                         STOP1165: !                        call abort
1144:                 endif1166:                 endif
1145:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i1167:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1146:         enddo1168:         enddo
1147: 1169: 
1148: 1170: 
1149:            tempN=01171:            tempN=0
1150: 1172: 
1151: 1173: 
1152: !$OMP PARALLEL1174: !$OMP PARALLEL
1153: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,C12,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  1175: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,C12,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  
1185:                  f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)1207:                  f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)
1186:                  if(r2 .gt. DHswitch2)then1208:                  if(r2 .gt. DHswitch2)then
1187:                   RD1=r1-DHswitch1209:                   RD1=r1-DHswitch
1188:                   f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))1210:                   f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))
1189:                   energy=energy+COEF2*RD1**2+COEF3*RD1**31211:                   energy=energy+COEF2*RD1**2+COEF3*RD1**3
1190:                  elseif(r2 .lt. DHswitch2)then1212:                  elseif(r2 .lt. DHswitch2)then
1191:                 ! normal DH term1213:                 ! normal DH term
1192:                  else1214:                  else
1193:                 ! things should have fallen in one of the previous two...1215:                 ! things should have fallen in one of the previous two...
1194:                   write(*,*) 'something went wrong with DH switching function'1216:                   write(*,*) 'something went wrong with DH switching function'
1195:                   STOP1217: !                call abort
1196:                  endif1218:                  endif
1197:                  f_over_r=f_over_r*sqrt(rm2)1219:                  f_over_r=f_over_r*sqrt(rm2)
1198:                  grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx1220:                  grad(C1*3-2) = grad(C1*3-2) + f_over_r * dx
1199:                  grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy1221:                  grad(C1*3-1) = grad(C1*3-1) + f_over_r * dy
1200:                  grad(C1*3)   = grad(C1*3)   + f_over_r * dz1222:                  grad(C1*3)   = grad(C1*3)   + f_over_r * dz
1201:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx1223:                  grad(C2*3-2) = grad(C2*3-2) - f_over_r * dx
1202:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy1224:                  grad(C2*3-1) = grad(C2*3-1) - f_over_r * dy
1203:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz1225:                  grad(C2*3)   = grad(C2*3)   - f_over_r * dz
1204: 1226: 
1205:                   endif1227:                   endif


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0