From 57a36a1183cb12334e9bc9ceacabbf038729d459 Mon Sep 17 00:00:00 2001 From: Adam Kazimierz Sieradzan Date: Thu, 24 Jul 2014 10:22:58 -0400 Subject: [PATCH] After long DEBUG of energy function --- source/unres/src_MD-M/elecont.f | 54 ++++++++++- source/unres/src_MD-M/energy_p_new-sep.F | 112 +++++++++++++++++++--- source/unres/src_MD-M/energy_p_new-sep_barrier.F | 17 +++- source/unres/src_MD-M/energy_p_new_barrier.F | 91 +++++++++++------- 4 files changed, 218 insertions(+), 56 deletions(-) diff --git a/source/unres/src_MD-M/elecont.f b/source/unres/src_MD-M/elecont.f index 4faf3e5..73325f2 100644 --- a/source/unres/src_MD-M/elecont.f +++ b/source/unres/src_MD-M/elecont.f @@ -52,6 +52,12 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ xmedi=xi+0.5*dxi ymedi=yi+0.5*dyi zmedi=zi+0.5*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize do 4 j=i+2,nct-1 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4 ind=ind+1 @@ -66,9 +72,49 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ dxj=c(1,j+1)-c(1,j) dyj=c(2,j+1)-c(2,j) dzj=c(3,j+1)-c(3,j) - xj=c(1,j)+0.5*dxj-xmedi - yj=c(2,j)+0.5*dyj-ymedi - zj=c(3,j)+0.5*dzj-zmedi + xj=c(1,j)+0.5*dxj + yj=c(2,j)+0.5*dyj + zj=c(3,j)+0.5*dzj + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif + rij=xj*xj+yj*yj+zj*zj + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) rrmij=1.0/(xj*xj+yj*yj+zj*zj) rmij=sqrt(rrmij) r3ij=rrmij*rmij @@ -94,7 +140,7 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ econt(ncont)=eesij endif ees=ees+eesij - evdw=evdw+evdwij + evdw=evdw+evdwij*sss 4 continue 1 continue if (lprint) then diff --git a/source/unres/src_MD-M/energy_p_new-sep.F b/source/unres/src_MD-M/energy_p_new-sep.F index 85f3e1f..de386c6 100644 --- a/source/unres/src_MD-M/energy_p_new-sep.F +++ b/source/unres/src_MD-M/energy_p_new-sep.F @@ -2128,6 +2128,12 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxysize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -2143,13 +2149,50 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi - rij=xj*xj+yj*yj+zj*zj + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif + rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) sss=sscale(rij/rpp(iteli,itelj)) + sssgrad=sscagrad(sqrt(rij)) if (sss.lt.1.0d0) then rmij=1.0D0/rij r3ij=rrmij*rmij @@ -2168,9 +2211,9 @@ C C Calculate contributions to the Cartesian gradient. C facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss) - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj do k=1,3 ghalf=0.5D0*ggg(k) @@ -2227,6 +2270,12 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxysize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -2242,13 +2291,50 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) sss=sscale(rij/rpp(iteli,itelj)) + sssgrad=sscagrad(sqrt(rij)) if (sss.gt.0.0d0) then rmij=1.0D0/rij r3ij=rrmij*rmij @@ -2267,9 +2353,9 @@ C C Calculate contributions to the Cartesian gradient. C facvdw=-6*rrmij*(ev1+evdwij)*sss - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj do k=1,3 ghalf=0.5D0*ggg(k) diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index a080b08..7d2d27f 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -1311,6 +1311,7 @@ cd write(iout,*) 'EE',EE(:,:,i) cd enddo cd call check_vecgrad cd stop +C print *,"WCHODZE3" if (icheckgrad.eq.1) then do i=1,nres-1 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i))) @@ -1406,6 +1407,12 @@ C xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=num_cont_hb(i) call eelecij_scale(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) @@ -1496,6 +1503,7 @@ C 13-go grudnia roku pamietnego... & 0.0d0,0.0d0,1.0d0/ c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j +C print *,"WCHODZE2" ind=ind+1 iteli=itel(i) itelj=itel(j) @@ -1519,7 +1527,7 @@ cd write (iout,*) "eelecij",i,j if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj zj_safe=zj @@ -1530,7 +1538,7 @@ cd write (iout,*) "eelecij",i,j xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj @@ -2164,6 +2172,7 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions double precision scal_el /0.5d0/ #endif evdw1=0.0D0 +C print *,"WCHODZE" c write (iout,*) "iatel_s_vdw",iatel_s_vdw, c & " iatel_e_vdw",iatel_e_vdw call flush(iout) @@ -2211,7 +2220,7 @@ c & ' ielend',ielend_vdw(i) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj zj_safe=zj @@ -2222,7 +2231,7 @@ c & ' ielend',ielend_vdw(i) xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 23aafae..ad328c0 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2010,7 +2010,7 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' dimension ggg(3) -cd write(iout,*) 'In EELEC_soft_sphere' +C write(iout,*) 'In EELEC_soft_sphere' ees=0.0D0 evdw1=0.0D0 eel_loc=0.0d0 @@ -2025,6 +2025,12 @@ cd write(iout,*) 'In EELEC_soft_sphere' xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -2038,10 +2044,49 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) dxj=dc(1,j) dyj=dc(2,j) dzj=dc(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif rij=xj*xj+yj*yj+zj*zj + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) if (rij.lt.r0ijsq) then evdw1ij=0.25d0*(rij-r0ijsq)**2 fac=rij-r0ijsq @@ -2049,13 +2094,13 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) evdw1ij=0.0d0 fac=0.0d0 endif - evdw1=evdw1+evdw1ij + evdw1=evdw1+evdw1ij*sss C C Calculate contributions to the Cartesian gradient. C - ggg(1)=fac*xj - ggg(2)=fac*yj - ggg(3)=fac*zj + ggg(1)=fac*xj*sssgrad + ggg(2)=fac*yj*sssgrad + ggg(3)=fac*zj*sssgrad do k=1,3 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) @@ -2928,31 +2973,6 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi -C Return atom into box, boxxsize is size of box in x dimension -c 184 continue -c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize -c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize -C Condition for being inside the proper box -c if ((xmedi.gt.((0.5d0)*boxxsize)).or. -c & (xmedi.lt.((-0.5d0)*boxxsize))) then -c go to 184 -c endif -c 185 continue -c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize -c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize -cC Condition for being inside the proper box -c if ((ymedi.gt.((0.5d0)*boxysize)).or. -c & (ymedi.lt.((-0.5d0)*boxysize))) then -c go to 185 -c endif -c 186 continue -c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize -c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize -cC Condition for being inside the proper box -c if ((zmedi.gt.((0.5d0)*boxzsize)).or. -c & (zmedi.lt.((-0.5d0)*boxzsize))) then -c go to 186 -c endif xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) @@ -3169,7 +3189,8 @@ C zj=c(3,j)+0.5D0*dzj-zmedi if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ" + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj zj_safe=zj @@ -3180,7 +3201,7 @@ C zj=c(3,j)+0.5D0*dzj-zmedi xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj -- 1.7.9.5