From 46f3f544d5a60812e2e11e03b619903f2ed052eb Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 29 Mar 2017 13:17:46 +0200 Subject: [PATCH] small chanegs in nanotube + working wham for lipid --- source/unres/src_MD-M/COMMON.TORSION | 7 ++-- source/unres/src_MD-M/energy_p_new_barrier.F | 50 +++++++++++++------------- source/wham/src-M/COMMON.SHIELD | 6 ++-- source/wham/src-M/energy_p_new.F | 40 ++++++++++++++++++--- source/wham/src-M/parmread.F | 1 + 5 files changed, 70 insertions(+), 34 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION index e7c54a0..20cf153 100644 --- a/source/unres/src_MD-M/COMMON.TORSION +++ b/source/unres/src_MD-M/COMMON.TORSION @@ -40,14 +40,15 @@ C surface double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde, & b,bnew1,bnew2,eenew,gtb1,gtb2,eeold,gtee integer nloctyp,iloctyp(-ntyp1:ntyp1),itype2loc(-ntyp1:ntyp1) - common/fourier/ b1(2,maxres),b2(2,maxres),b(13,-ntyp:ntyp), + common/fourier/ b1(2,-maxres:maxres),b2(2,-maxres:maxres), + & b(13,-ntyp:ntyp), & bnew1(3,2,-ntyp:ntyp),bnew2(3,2,-ntyp:ntyp), & cc(2,2,-ntyp:ntyp), & dd(2,2,-ntyp:ntyp),eeold(2,2,-ntyp:ntyp), & eenew(2,-ntyp:ntyp), & ee(2,2,maxres), & ctilde(2,2,-ntyp:ntyp), - & dtilde(2,2,-ntyp:ntyp),b1tilde(2,maxres), - & b2tilde(2,maxres), + & dtilde(2,2,-ntyp:ntyp),b1tilde(2,-maxres:maxres), + & b2tilde(2,-maxres:maxres), & gtb1(2,maxres),gtb2(2,maxres),gtEE(2,2,maxres), & nloctyp,iloctyp,itype2loc 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 1390fe3..e3a17dd 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -4647,8 +4647,8 @@ c & a33*gmuji2(4) #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'eelloc',i,j,eel_loc_ij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2f7.3)') + & 'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2) c if (eel_loc_ij.ne.0) c & write (iout,'(a4,2i4,8f9.5)')'chuj', c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) @@ -12611,17 +12611,17 @@ C now calculate distance from center of tube and direction vectors zmin=boxzsize do j=-1,1 - vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) vectube(1)=vectube(1)+boxxsize*j - vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize) vectube(2)=vectube(2)+boxysize*j - vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize) + vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize) vectube(3)=vectube(3)+boxzsize*j - xminact=abs(vectube(1)-tubecenter(1)) - yminact=abs(vectube(2)-tubecenter(2)) - zminact=abs(vectube(3)-tubecenter(3)) + xminact=dabs(vectube(1)-tubecenter(1)) + yminact=dabs(vectube(2)-tubecenter(2)) + zminact=dabs(vectube(3)-tubecenter(3)) if (xmin.gt.xminact) then xmin=xminact @@ -12674,13 +12674,13 @@ C go to 667 enecavtube(i)=0.0 faccav=0.0 else - denominator=(1.0+dcavtubpep*rdiff6*rdiff6) + denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6) enecavtube(i)= - & (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep) + & (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep) & /denominator enecavtube(i)=0.0 - faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff)) - & *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(rdiff) + faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff)) + & *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff) & +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0) & /denominator**2.0d0 C faccav=0.0 @@ -12701,7 +12701,7 @@ C now direction of gg_tube vector enddo do i=itube_start,itube_end - enecavtube(i)=0.0 + enecavtube(i)=0.0d0 C Lets not jump over memory as we use many times iti iti=itype(i) C lets ommit dummy atoms for now @@ -12713,17 +12713,17 @@ C .or.(iti.eq.10) ymin=boxysize zmin=boxzsize do j=-1,1 - vectube(1)=mod((c(1,i+nres)),boxxsize) + vectube(1)=dmod((c(1,i+nres)),boxxsize) vectube(1)=vectube(1)+boxxsize*j - vectube(2)=mod((c(2,i+nres)),boxysize) + vectube(2)=dmod((c(2,i+nres)),boxysize) vectube(2)=vectube(2)+boxysize*j - vectube(3)=mod((c(3,i+nres)),boxzsize) + vectube(3)=dmod((c(3,i+nres)),boxzsize) vectube(3)=vectube(3)+boxzsize*j - xminact=abs(vectube(1)-tubecenter(1)) - yminact=abs(vectube(2)-tubecenter(2)) - zminact=abs(vectube(3)-tubecenter(3)) + xminact=dabs(vectube(1)-tubecenter(1)) + yminact=dabs(vectube(2)-tubecenter(2)) + zminact=dabs(vectube(3)-tubecenter(3)) if (xmin.gt.xminact) then xmin=xminact @@ -12772,16 +12772,16 @@ C now direction of gg_tube vector C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12) if (acavtub(iti).eq.0.0d0) then C go to 667 - enecavtube(i+nres)=0.0 - faccav=0.0 + enecavtube(i+nres)=0.0d0 + faccav=0.0d0 else - denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6) + denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6) enecavtube(i+nres)= - & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti)) + & (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti)) & /denominator C enecavtube(i)=0.0 - faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff)) - & *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff) + faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff)) + & *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff) & +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0) & /denominator**2.0d0 C faccav=0.0 diff --git a/source/wham/src-M/COMMON.SHIELD b/source/wham/src-M/COMMON.SHIELD index 1f96c94..ac7ac98 100644 --- a/source/wham/src-M/COMMON.SHIELD +++ b/source/wham/src-M/COMMON.SHIELD @@ -1,9 +1,11 @@ double precision VSolvSphere,VSolvSphere_div,long_r_sidechain, & short_r_sidechain,fac_shield,grad_shield_side,grad_shield, - & buff_shield,wshield,grad_shield_loc + & buff_shield,wshield,grad_shield_loc,lipscale,sslipi,sslipj, + & ssgradlipi,ssgradlipj integer ishield_list,shield_list,ees0plist common /shield/ VSolvSphere,VSolvSphere_div,buff_shield, - & long_r_sidechain(ntyp), + & long_r_sidechain(ntyp),sslipi,sslipj, + & ssgradlipi,ssgradlipj,lipscale, & short_r_sidechain(ntyp),fac_shield(maxres),wshield, & grad_shield_side(3,maxcont,-1:maxres),grad_shield(3,-1:maxres), & grad_shield_loc(3,maxcont,-1:maxres), diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 378eaae..a36e177 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -2398,8 +2398,12 @@ C fac_shield(j)=0.6 fac_shield(j)=1.0 eesij=(el1+el2) ees=ees+eesij + &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + endif evdw1=evdw1+evdwij*sss + &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') c &'evdw1',i,j,evdwij c &,iteli,itelj,aaa,evdw1 @@ -2414,7 +2418,11 @@ C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE facvdw=-6*rrmij*(ev1+evdwij)*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facel=-3*rrmij*(el1+eesij) + &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + fac1=fac erij(1)=xj*rmij erij(2)=yj*rmij @@ -2498,8 +2506,13 @@ C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj if (sss.gt.0.0) then ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + else ggg(1)=0.0 ggg(2)=0.0 @@ -2510,6 +2523,11 @@ C ggg(3)=facvdw*zj gvdwpp(k,i)=gvdwpp(k,i)+ghalf gvdwpp(k,j)=gvdwpp(k,j)+ghalf enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 + * * Loop over residues i+1 thru j-1. * @@ -2520,6 +2538,8 @@ C ggg(3)=facvdw*zj enddo #else facvdw=(ev1+evdwij)*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facel=el1+eesij fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) @@ -2564,6 +2584,8 @@ cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) & *fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + enddo do k=1,3 ghalf=0.5D0*ggg(k) @@ -2571,11 +2593,15 @@ cd & (dcosg(k),k=1,3) & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) & *fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + gelc(k,j)=gelc(k,j)+ghalf & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) & *fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + enddo do k=i+1,j-1 do l=1,3 @@ -2831,9 +2857,11 @@ C fac_shield(j)=0.6 eel_loc_ij=eel_loc_ij & *fac_shield(i)*fac_shield(j) &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) -c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij -C write (iout,'(a6,2i5,0pf7.3)') -C & 'eelloc',i,j,eel_loc_ij +C write (iout,'(a3,i4,a3,i4,a8,4f8.3)') +C & 'i',i,' j',j,' eel_loc_ij',eel_loc_ij,sslipi, +C & sslipj,lipscale +C write (iout,'(a6,2i5,0pf7.3,2f7.3)') +C & 'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2) C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) C eel_loc=eel_loc+eel_loc_ij @@ -3272,7 +3300,9 @@ C fac_shield(j)=0.6 eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) - + write (iout,'(a3,i4,a3,i4,a8,4f8.3)') + & 'i',i,' j',j,' eturn3',eello_t3,sslipi, + & sslipj,lipscale cd write (2,*) 'i,',i,' j',j,'eello_turn3', cd & 0.5d0*(pizda(1,1)+pizda(2,2)), cd & ' eello_turn3_num',4*eello_turn3_num @@ -3431,6 +3461,8 @@ C fac_shield(j)=0.6 eello_turn4=eello_turn4-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index 2fc37d3..ce1bf31 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -71,6 +71,7 @@ c If reading not own parameters, skip assignment call reada(controlcard,"ATRISS",atriss,0.3D0) call reada(controlcard,"BTRISS",btriss,0.02D0) call reada(controlcard,"CTRISS",ctriss,1.0D0) + call reada(controlcard,"LIPSCALE",lipscale,1.3D0) dyn_ss=(index(controlcard,'DYN_SS').gt.0) write(iout,*) "ATRISS",atriss write(iout,*) "BTRISS",btriss -- 1.7.9.5