From 8c90e4c2c3004fdd391d6540f654ef9a1b9b49aa Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 12 May 2014 19:15:54 +0200 Subject: [PATCH] po przyspieszeniu przed checkgradem i sympcheckiem --- source/unres/src_MD-M/energy_p_new_barrier.F | 235 ++++++++++++++++++++------ 1 file changed, 181 insertions(+), 54 deletions(-) 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 252cff0..c8acad2 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -152,9 +152,9 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" eello_turn4=0.0d0 endif else -c write (iout,*) "Soft-spheer ELEC potential" - call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, - & eello_turn4) + write (iout,*) "Soft-spheer ELEC potential" +c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, +c & eello_turn4) endif c print *,"Processor",myrank," computed UELEC" C @@ -1425,9 +1425,9 @@ c if (icall.eq.0) lprn=.false. ind=0 C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0 C we have the original box) - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle @@ -1466,9 +1466,9 @@ c endif if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - xi=xi+xshift*boxxsize - yi=yi+yshift*boxysize - zi=zi+zshift*boxzsize +C xi=xi+xshift*boxxsize +C yi=yi+yshift*boxysize +C zi=zi+zshift*boxzsize dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) @@ -1544,12 +1544,43 @@ c endif 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 + subchap=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 + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) - xj=xj-xi - yj=yj-yi - zj=zj-zi +C xj=xj-xi +C yj=yj-yi +C zj=zj-zi c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi c write (iout,*) "j",j," dc_norm", c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) @@ -1623,9 +1654,9 @@ C Calculate angular part of the gradient. enddo ! j enddo ! iint enddo ! i - enddo ! zshift - enddo ! yshift - enddo ! xshift +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return @@ -2989,9 +3020,9 @@ c endif num_cont_hb(i)=num_conti enddo ! i C Loop over all neighbouring boxes - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c @@ -3015,9 +3046,9 @@ c if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize - xmedi=xmedi+xshift*boxxsize - ymedi=ymedi+yshift*boxysize - zmedi=zmedi+zshift*boxzsize +C xmedi=xmedi+xshift*boxxsize +C ymedi=ymedi+yshift*boxysize +C zmedi=zmedi+zshift*boxzsize C Return tom into box, boxxsize is size of box in x dimension c 164 continue @@ -3057,9 +3088,9 @@ c write (iout,*) i,j,itype(i),itype(j) enddo ! j num_cont_hb(i)=num_conti enddo ! i - enddo ! zshift - enddo ! yshift - enddo ! xshift +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift c write (iout,*) "Number of loop steps in EELEC:",ind cd do i=1,nres @@ -3138,7 +3169,37 @@ 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 + 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 C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC c 174 continue c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize @@ -3165,9 +3226,9 @@ c & (zj.lt.((-0.5d0)*boxzsize))) then c go to 176 c endif C endif !endPBC condintion - xj=xj-xmedi - yj=yj-ymedi - zj=zj-zmedi +C xj=xj-xmedi +C yj=yj-ymedi +C zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj sss=sscale(sqrt(rij)) @@ -4110,9 +4171,9 @@ C r0_scp=4.5d0 cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) @@ -4150,9 +4211,9 @@ c endif if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - xi=xi+xshift*boxxsize - yi=yi+yshift*boxysize - zi=zi+zshift*boxzsize +C xi=xi+xshift*boxxsize +C yi=yi+yshift*boxysize +C zi=zi+zshift*boxzsize do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -4195,10 +4256,41 @@ c go to 176 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 + subchap=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 + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif c c endif - xj=xj-xi - yj=yj-yi - zj=zj-zi +C xj=xj-xi +C yj=yj-yi +C zj=zj-zi rij=xj*xj+yj*yj+zj*zj r0ij=r0_scp @@ -4251,9 +4343,9 @@ cgrad enddo enddo ! iint enddo ! i - enddo !zshift - enddo !yshift - enddo !xshift +C enddo !zshift +C enddo !yshift +C enddo !xshift return end C----------------------------------------------------------------------------- @@ -4278,11 +4370,12 @@ C dimension ggg(3) evdw2=0.0D0 evdw2_14=0.0d0 +c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) @@ -4295,9 +4388,10 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - xi=xi+xshift*boxxsize - yi=yi+yshift*boxysize - zi=zi+zshift*boxzsize +c xi=xi+xshift*boxxsize +c yi=yi+yshift*boxysize +c zi=zi+zshift*boxzsize +c print *,xi,yi,zi,'polozenie i' C Return atom into box, boxxsize is size of box in x dimension c 134 continue c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize @@ -4368,13 +4462,46 @@ c if ((zj.gt.((0.5d0)*boxzsize)).or. c & (zj.lt.((-0.5d0)*boxzsize))) then c go to 176 c endif - xj=xj-xi - yj=yj-yi - zj=zj-zi +CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=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 + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif +c print *,xj,yj,zj,'polozenie j' rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +c print *,rrij sss=sscale(1.0d0/(dsqrt(rrij))) +c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz' +c if (sss.eq.0) print *,'czasem jest OK' + if (sss.le.0.0d0) cycle sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) - if (sss.gt.0.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) @@ -4427,14 +4554,14 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - endif !endif for sscale cutoff +c endif !endif for sscale cutoff enddo ! j enddo ! iint enddo ! i - enddo !zshift - enddo !yshift - enddo !xshift +c enddo !zshift +c enddo !yshift +c enddo !xshift do i=1,nct do j=1,3 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) -- 1.7.9.5