From 36865e3203724663beeb1f60abdbafbec07b1019 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Fri, 1 Apr 2016 09:10:54 +0200 Subject: [PATCH] AFM+tube+kcc --- PARAM/tube.parm | 50 +++--- source/unres/src_MD-M/COMMON.CHAIN | 6 +- source/unres/src_MD-M/COMMON.INTERACT | 4 +- source/unres/src_MD-M/energy_p_new_barrier.F | 210 +++++++++++++++++++++++++- source/unres/src_MD-M/parmread.F | 4 +- source/unres/src_MD-M/readrtns_CSA.F | 15 +- 6 files changed, 249 insertions(+), 40 deletions(-) diff --git a/PARAM/tube.parm b/PARAM/tube.parm index d1783aa..182c68f 100644 --- a/PARAM/tube.parm +++ b/PARAM/tube.parm @@ -1,26 +1,26 @@ - 200.00 2.20 - 1.00 0.50 - 1.00 0.20 - 1.00 0.30 - 1.00 1.00 - 1.00 1.00 - 200.00 2.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 - 1.00 1.00 + 200.00 2.20 3.0 + 1.00 0.50 3.0 + 1.00 0.20 3.0 + 1.00 0.30 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 200.00 2.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 + 1.00 1.00 3.0 diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index 8443ef4..372fa85 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -19,11 +19,13 @@ double precision boxxsize,boxysize,boxzsize,enecut,sscut, & sss,sssgrad, & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick, - & tubecenter,tubeR0 + & tubecenter,tubeR0, + & buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad, & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick common /afm/ forceAFMconst, distafminit,afmend,afmbeg, & velAFMconst, & totTafm - common /tube/ tubecenter(3),tubeR0 + common /tube/ tubecenter(3),tubeR0, + & buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT index c43ef6a..14d92ef 100644 --- a/source/unres/src_MD-M/COMMON.INTERACT +++ b/source/unres/src_MD-M/COMMON.INTERACT @@ -41,6 +41,8 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters. & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp), & distchainmax,nbondterm(ntyp) C 01/29/15 Lipidic parameters - double precision pepliptran,liptranene + double precision pepliptran,liptranene, + &tubetranene, tubetranenepep common /lipid/ pepliptran,liptranene(ntyp) + common /tubepar/ tubetranene(ntyp), tubetranenepep 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 e86bb6e..21c0197 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -307,9 +307,11 @@ C print *,"za lipidami" else if (selfguide.gt.0) then call AFMvel(Eafmforce) endif - if (TUBElog.gt.0) then + if (TUBElog.eq.1) then C print *,"just before call" call calctube(Etube) + elseif (TUBElog.eq.2) then + call calctube2(Etube) else Etube=0.0d0 endif @@ -11725,8 +11727,13 @@ C for UNRES C lets ommit dummy atoms for now if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle C now calculate distance from center of tube and direction vectors - vectube(1)=(c(1,i)+c(1,i+1))/2.0d0-tubecenter(1) - vectube(2)=(c(2,i)+c(2,i+1))/2.0d0-tubecenter(2) + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) @@ -11768,8 +11775,15 @@ C lets ommit dummy atoms for now C in UNRES uncomment the line below as GLY has no side-chain... C .or.(iti.eq.10) & ) cycle - vectube(1)=c(1,i+nres)-tubecenter(1) - vectube(2)=c(2,i+nres)-tubecenter(2) + vectube(1)=c(1,i+nres) + vectube(1)=mod(vectube(1),boxxsize) + if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize + vectube(2)=c(2,i+nres) + vectube(2)=mod(vectube(2),boxxsize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) C as the tube is infinity we do not calculate the Z-vector use of Z C as chosen axis @@ -11810,3 +11824,189 @@ C 4) add reading the center of tube C 5) add COMMONs C 6) add to zerograd +C----------------------------------------------------------------------- +C----------------------------------------------------------- +C This subroutine is to mimic the histone like structure but as well can be +C utilizet to nanostructures (infinit) small modification has to be used to +C make it finite (z gradient at the ends has to be changes as well as the x,y +C gradient has to be modified at the ends +C The energy function is Kihara potential +C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6) +C 4eps is depth of well sigma is r_minimum r is distance from center of tube +C and r0 is the excluded size of nanotube (can be set to 0 if we want just a +C simple Kihara potential + subroutine calctube2(Etube) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' + double precision tub_r,vectube(3),enetube(maxres*2) + Etube=0.0d0 + do i=1,2*nres + enetube(i)=0.0d0 + enddo +C first we calculate the distance from tube center +C first sugare-phosphate group for NARES this would be peptide group +C for UNRES + do i=1,nres +C lets ommit dummy atoms for now + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +C now calculate distance from center of tube and direction vectors + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + +C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) + +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6 +C write(iout,*) "TU13",i,rdiff6,enetube(i) +C print *,rdiff,rdiff6,pep_aa_tube +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6+ + & 6.0d0*pep_bb_tube)/rdiff6/rdiff +C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +C &rdiff,fac + +C now direction of gg_tube vector + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + enddo +C basically thats all code now we split for side-chains (REMEMBER to sum up at the END) + do i=1,nres +C Lets not jump over memory as we use many times iti + iti=itype(i) +C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) +C in UNRES uncomment the line below as GLY has no side-chain... + & .or.(iti.eq.10) + & ) cycle + vectube(1)=c(1,i+nres) + vectube(1)=mod(vectube(1),boxxsize) + if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize + vectube(2)=c(2,i+nres) + vectube(2)=mod(vectube(2),boxxsize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) +C THIS FRAGMENT MAKES TUBE FINITE + positi=(mod(c(3,i+nres),boxzsize)) + if (positi.le.0) positi=positi+boxzsize +C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop +c for each residue check if it is in lipid or lipid water border area +C respos=mod(c(3,i+nres),boxzsize) + print *,positi,bordtubebot,buftubebot,bordtubetop + if ((positi.gt.bordtubebot) + & .and.(positi.lt.bordtubetop)) then +C the energy transfer exist + if (positi.lt.buftubebot) then + fracinbuf=1.0d0- + & ((positi-bordtubebot)/tubebufthick) +C lipbufthick is thickenes of lipid buffore + sstube=sscalelip(fracinbuf) + ssgradtube=-sscagradlip(fracinbuf)/tubebufthick + print *,ssgradtube, sstube,tubetranene(itype(i)) + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) + gg_tube_SC(3,i)=gg_tube_SC(3,i) + &+ssgradtube*tubetranene(itype(i)) + gg_tube(3,i-1)= gg_tube(3,i-1) + &+ssgradtube*tubetranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.buftubetop) then + fracinbuf=1.0d0- + &((bordtubetop-positi)/tubebufthick) + sstube=sscalelip(fracinbuf) + ssgradtube=sscagradlip(fracinbuf)/tubebufthick + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + sstube=1.0d0 + ssgradtube=0.0d0 + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C print *,"I am in true lipid" + endif + else +C sstube=0.0d0 +C ssgradtube=0.0d0 + cycle + endif ! if in lipid or buffor +CEND OF FINITE FRAGMENT +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6) + & *sstube+enetube(i+nres) +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+ + & 6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube +C now direction of gg_tube vector + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + gg_tube_SC(3,i)=gg_tube_SC(3,i) + &+ssgradtube*enetube(i+nres)/sstube + gg_tube(3,i-1)= gg_tube(3,i-1) + &+ssgradtube*enetube(i+nres)/sstube + + enddo + do i=1,2*nres + Etube=Etube+enetube(i) + enddo +C print *,"ETUBE", etube + return + end +C TO DO 1) add to total energy +C 2) add to gradient summation +C 3) add reading parameters (AND of course oppening of PARAM file) +C 4) add reading the center of tube +C 5) add COMMONs +C 6) add to zerograd + diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index a38147f..5001b6b 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -1367,7 +1367,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) epspeptube=dabs(epspeptube) pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube - write(iout,*) pep_aa_tube,pep_bb_tube + write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep do i=1,ntyp read(itube,*) epssctube,sigmasctube sigmasctube=sigmasctube**6 @@ -1375,7 +1375,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) epssctube=dabs(epssctube) sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube - write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i) + write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i) enddo #ifdef OLDSCP diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 89ddeac..e40717f 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -147,11 +147,6 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'TUBEMOD',tubelog,0) write (iout,*) TUBElog,"TUBEMODE" - if (TUBElog.gt.0) then - call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) - call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) - call reada(controlcard,"RTUBE",tubeR0,0.0d0) - endif call readi(controlcard,'IPRINT',iprint,0) C SHIELD keyword sets if the shielding effect of side-chains is used C 0 denots no shielding is used all peptide are equally despite the @@ -270,6 +265,16 @@ C endif write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot write (iout,*) "SHIELD MODE",shield_mode + if (TUBElog.gt.0) then + call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) + call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) + call reada(controlcard,"RTUBE",tubeR0,0.0d0) + call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) + call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) + call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) + buftubebot=bordtubebot+tubebufthick + buftubetop=bordtubetop-tubebufthick + endif if (shield_mode.gt.0) then pi=3.141592d0 C VSolvSphere the volume of solving sphere -- 1.7.9.5