hdiff output

r30596/SBM.f 2016-07-06 15:35:20.011014640 +0100 r30595/SBM.f 2016-07-06 15:35:20.379019361 +0100
257: ! read in position restraints257: ! read in position restraints
258:        read(30,*) 258:        read(30,*) 
259:        read(30,*) SBMPRN259:        read(30,*) SBMPRN
260:        write(*,*) SBMPRN, ' position restraints'260:        write(*,*) SBMPRN, ' position restraints'
261:        do I=1,SBMPRN261:        do I=1,SBMPRN
262:             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),
263:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6),263:      Q         SBMPRK(I,4),SBMPRK(I,5),SBMPRK(I,6),
264:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)264:      Q         SBMPRX(I,1),SBMPRX(I,2),SBMPRX(I,3)
265:             if(TARR(SBMPRI(I)) .ne. 0)then265:             if(TARR(SBMPRI(I)) .ne. 0)then
266:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)266:                 write(*,*) 'more than one restraint provided for atom ',SBMPRI(I)
267: !                call abort267:                 call abort
268:             endif268:             endif
269:             TARR(SBMPRI(I))=1269:             TARR(SBMPRI(I))=1
270:         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  
271:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'271:           write(*,*) 'FATAL ERROR: Non-zero restraint cross-terms not supported'
272:         endif272:         endif
273:        enddo 273:        enddo 
274: 274: 
275:        close(30)275:        close(30)
276:        end276:        end
277: 277: 
900:         enddo900:         enddo
901: 901: 
902:         maxgridX=int((maxX-minX)/gridsize)+1902:         maxgridX=int((maxX-minX)/gridsize)+1
903:         maxgridY=int((maxY-minY)/gridsize)+1903:         maxgridY=int((maxY-minY)/gridsize)+1
904:         maxgridZ=int((maxZ-minZ)/gridsize)+1904:         maxgridZ=int((maxZ-minZ)/gridsize)+1
905: 905: 
906:         if(maxgridX .ge. maxgrid .or. 906:         if(maxgridX .ge. maxgrid .or. 
907:      Q  maxgridY .ge. maxgrid .or.907:      Q  maxgridY .ge. maxgrid .or.
908:      Q  maxgridZ .ge. maxgrid )then908:      Q  maxgridZ .ge. maxgrid )then
909:         write(*,*) 'system got too big for grid searching...'909:         write(*,*) 'system got too big for grid searching...'
910: !        call abort910:         call abort
911:         endif911:         endif
912: 912: 
913: 913: 
914: !!  Add a second grid that only includes atoms that are charged914: !!  Add a second grid that only includes atoms that are charged
915: 915: 
916:         do i=1,maxgrid916:         do i=1,maxgrid
917:          do j=1,maxgrid917:          do j=1,maxgrid
918:           do k=1,maxgrid918:           do k=1,maxgrid
919:                 Ngrid(i,j,k)=0919:                 Ngrid(i,j,k)=0
920:           enddo920:           enddo
922:         enddo922:         enddo
923:         do i=1,NATOMS923:         do i=1,NATOMS
924: 924: 
925:                 Xgrid=int((X(i)-minX)/gridsize)+1925:                 Xgrid=int((X(i)-minX)/gridsize)+1
926:                 Ygrid=int((Y(i)-minY)/gridsize)+1926:                 Ygrid=int((Y(i)-minY)/gridsize)+1
927:                 Zgrid=int((Z(i)-minZ)/gridsize)+1927:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
928:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1928:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
929:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then929:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
930:                         write(*,*) 'too many atoms in a grid'930:                         write(*,*) 'too many atoms in a grid'
931:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid931:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
932: !                        call abort932:                         call abort
933:                 endif933:                 endif
934:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i934:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
935:         enddo935:         enddo
936: 936: 
937: 937: 
938:            tempN=0938:            tempN=0
939: 939: 
940: ! add a second loop that goes over charged atoms940: ! add a second loop that goes over charged atoms
941: 941: 
942: !$OMP PARALLEL942: !$OMP PARALLEL
977:                 f_over_r=-STT*12.0*rm14977:                 f_over_r=-STT*12.0*rm14
978:                 RD1=sqrt(r2)-NCswitch978:                 RD1=sqrt(r2)-NCswitch
979:                 if(r2 .gt. Rswitch2)then979:                 if(r2 .gt. Rswitch2)then
980:                         f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)980:                         f_over_r=f_over_r+(SA*RD1**2+SB*RD1**3)*sqrt(rm2)
981:                         etemp=etemp+SA/3.0*RD1**3+SB/4.0*RD1**4981:                         etemp=etemp+SA/3.0*RD1**3+SB/4.0*RD1**4
982:                 elseif(r2 .lt. Rswitch2)then982:                 elseif(r2 .lt. Rswitch2)then
983:                 ! normal repulsive term983:                 ! normal repulsive term
984:                 else984:                 else
985:                 ! things should have fallen in one of the previous two...985:                 ! things should have fallen in one of the previous two...
986:                 write(*,*) 'something went wrong with switching function'986:                 write(*,*) 'something went wrong with switching function'
987: !                call abort987:                 call abort
988:                 endif988:                 endif
989:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx989:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx
990:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy990:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy
991:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz991:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz
992: 992: 
993:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx993:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx
994:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy994:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy
995:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz995:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz
996: 996: 
997:             endif997:             endif
1129:         enddo1129:         enddo
1130: 1130: 
1131:         maxgridX=int((maxX-minX)/gridsize)+11131:         maxgridX=int((maxX-minX)/gridsize)+1
1132:         maxgridY=int((maxY-minY)/gridsize)+11132:         maxgridY=int((maxY-minY)/gridsize)+1
1133:         maxgridZ=int((maxZ-minZ)/gridsize)+11133:         maxgridZ=int((maxZ-minZ)/gridsize)+1
1134: 1134: 
1135:         if(maxgridX .ge. maxgrid .or. 1135:         if(maxgridX .ge. maxgrid .or. 
1136:      Q  maxgridY .ge. maxgrid .or.1136:      Q  maxgridY .ge. maxgrid .or.
1137:      Q  maxgridZ .ge. maxgrid )then1137:      Q  maxgridZ .ge. maxgrid )then
1138:         write(*,*) 'system got too big for grid searching...'1138:         write(*,*) 'system got too big for grid searching...'
1139: !        call abort1139:         call abort
1140:         endif1140:         endif
1141: 1141: 
1142: 1142: 
1143: !!  Add a second grid that only includes atoms that are charged1143: !!  Add a second grid that only includes atoms that are charged
1144: 1144: 
1145:         do i=1,maxgrid1145:         do i=1,maxgrid
1146:          do j=1,maxgrid1146:          do j=1,maxgrid
1147:           do k=1,maxgrid1147:           do k=1,maxgrid
1148:                 Ngrid(i,j,k)=01148:                 Ngrid(i,j,k)=0
1149:           enddo1149:           enddo
1154:         do ii=2,SBMCHARGEON(1)1154:         do ii=2,SBMCHARGEON(1)
1155:                 i=SBMCHARGEON(ii)1155:                 i=SBMCHARGEON(ii)
1156: 1156: 
1157:                 Xgrid=int((X(i)-minX)/gridsize)+11157:                 Xgrid=int((X(i)-minX)/gridsize)+1
1158:                 Ygrid=int((Y(i)-minY)/gridsize)+11158:                 Ygrid=int((Y(i)-minY)/gridsize)+1
1159:                 Zgrid=int((Z(i)-minZ)/gridsize)+11159:                 Zgrid=int((Z(i)-minZ)/gridsize)+1
1160:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+11160:                 Ngrid(Xgrid,Ygrid,Zgrid)=Ngrid(Xgrid,Ygrid,Zgrid)+1
1161:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then1161:                 if(Ngrid(Xgrid,Ygrid,Zgrid) .gt. maxpergrid)then
1162:                         write(*,*) 'too many atoms in a grid'1162:                         write(*,*) 'too many atoms in a grid'
1163:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid1163:                         write(*,*) Ngrid(Xgrid,Ygrid,Zgrid),Xgrid,Ygrid,Zgrid
1164: !                        call abort1164:                         call abort
1165:                 endif1165:                 endif
1166:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i1166:                 grid(Xgrid,Ygrid,Zgrid,Ngrid(Xgrid,Ygrid,Zgrid))=i
1167:         enddo1167:         enddo
1168: 1168: 
1169: 1169: 
1170:            tempN=01170:            tempN=0
1171: 1171: 
1172: 1172: 
1173: !$OMP PARALLEL1173: !$OMP PARALLEL
1174: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  1174: !$OMP& PRIVATE(i,ii,iii,jj,kk,jjj,C1,C2,r2,rm2,rm14,f_over_r,RD1,dx,dy,dz,Xgrid,Ygrid,Zgrid,TID,C2T)  
1213:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)1213:                 f_over_r=PREFACTOR*C2T*SBMDHVP(kappa,r1)
1214:                 if(r2 .gt. DHswitch2)then1214:                 if(r2 .gt. DHswitch2)then
1215:                 RD1=r1-DHswitch1215:                 RD1=r1-DHswitch
1216:                         f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))1216:                         f_over_r=(f_over_r+C2T*(2*COEF2*RD1+3*COEF3*RD1**2))
1217:                         etemp=etemp+COEF2*RD1**2+COEF3*RD1**31217:                         etemp=etemp+COEF2*RD1**2+COEF3*RD1**3
1218:                 elseif(r2 .lt. DHswitch2)then1218:                 elseif(r2 .lt. DHswitch2)then
1219:                 ! normal DH term1219:                 ! normal DH term
1220:                 else1220:                 else
1221:                 ! things should have fallen in one of the previous two...1221:                 ! things should have fallen in one of the previous two...
1222:                 write(*,*) 'something went wrong with DH switching function'1222:                 write(*,*) 'something went wrong with DH switching function'
1223: !                call abort1223:                 call abort
1224:                 endif1224:                 endif
1225:                 f_over_r=f_over_r*sqrt(rm2)1225:                 f_over_r=f_over_r*sqrt(rm2)
1226:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx1226:               gradOMP(C1*3-2) = gradOMP(C1*3-2) + f_over_r * dx
1227:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy1227:               gradOMP(C1*3-1) = gradOMP(C1*3-1) + f_over_r * dy
1228:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz1228:               gradOMP(C1*3)   = gradOMP(C1*3)   + f_over_r * dz
1229: 1229: 
1230:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx1230:               gradOMP(C2*3-2) =  gradOMP(C2*3-2) - f_over_r * dx
1231:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy1231:               gradOMP(C2*3-1) =  gradOMP(C2*3-1) - f_over_r * dy
1232:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz1232:               gradOMP(C2*3)   =  gradOMP(C2*3)   - f_over_r * dz
1233: 1233: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0