--- /dev/null
+ 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
+
common /from_zscore/ nz_start,nz_end,iz_sc
double precision boxxsize,boxysize,boxzsize,enecut,sscut,
& sss,sssgrad,
- & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+ & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick,
+ & tubecenter,tubeR0
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
integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
& inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
- & shield_mode,tor_mode
+ & shield_mode,tor_mode,tubelog
logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
& sideadd,lsecondary,read_cart,unres_pdb,
& vdisulf,searchsc,lmuca,dccart,extconf,out1file,
& icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint,
& overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
& ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
- & selfguide,AFMlog,shield_mode,tor_mode,
+ & selfguide,AFMlog,shield_mode,tor_mode,tubelog,
& constr_dist,gnorm_check,gradout,split_ene,with_theta_constr,
& symetr
C... minim = .true. means DO minimization.
& gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gloc_x,dtheta,dphi,dalpha,
& domega,gscloc,gsclocx,gradcorr,gradcorr_long,gradcorr5_long,
& gradcorr6_long,gcorr6_turn_long,gvdwx,gshieldx,gradafm,
+ & gg_tube,gg_tube_SC,
& gshieldc, gshieldc_loc, gshieldx_ec, gshieldc_ec,
& gshieldc_loc_ec, gshieldx_t3,gshieldc_t3,gshieldc_loc_t3,
& gshieldx_t4, gshieldc_t4,gshieldc_loc_t4,gshieldx_ll,
& gshieldc_loc_t4(3,-1:maxres),
& gshieldx_ll(3,-1:maxres), gshieldc_ll(3,-1:maxres),
& gshieldc_loc_ll(3,-1:maxres),
- & gradafm(3,-1:maxres),
+ & gradafm(3,-1:maxres),gg_tube(3,-1:maxres),
+ & gg_tube_sc(3,-1:maxres),
& gradx_scp(3,-1:maxres),gvdwc_scp(3,-1:maxres),
& ghpbx(3,-1:maxres),
& ghpbc(3,-1:maxres),gloc(maxvar,2),gradcorr(3,-1:maxres),
C-----------------------------------------------------------------------
integer n_ene_comp,rescale_mode
common /ffield/ wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,
+ & wtube,
& wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
& wturn6,wvdwpp,weights(n_ene),wliptran,temp0,
& scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp,
double precision aa,bb,augm,aad,bad,app,bpp,ale6,ael3,ael6,
- &aa_lip,bb_lip,aa_aq,bb_aq
+ &aa_lip,bb_lip,aa_aq,bb_aq,sc_aa_tube_par,sc_bb_tube_par,
+ & pep_aa_tube,pep_bb_tube
+
integer expon,expon2
integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro,
& ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr,iscpstart,
& ispp,iscp
common /interact/aa_aq(ntyp,ntyp),bb_aq(ntyp,ntyp),
& aa_lip(ntyp,ntyp),bb_lip(ntyp,ntyp),
+ & sc_aa_tube_par(ntyp),sc_bb_tube_par(ntyp),
+ & pep_aa_tube,pep_bb_tube,
& augm(ntyp,ntyp),
& aad(ntyp,2),bad(ntyp,2),app(2,2),bpp(2,2),ael6(2,2),ael3(2,2),
& expon,expon2,nnt,nct,nint_gr(maxres),istart(maxres,maxint_gr),
& itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,istat,
& ientin,ientout,izs1,isecpred,ibond,irest2,iifrag,icart,
& irest1,isccor,ithep_pdb,irotam_pdb,itorkcc,ithetkcc,
- & iliptranpar
+ & iliptranpar,itube
common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,
& irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,
& istat,ientin,ientout,izs1,isecpred,ibond,irest2,iifrag,
& icart,irest1,isccor,ithep_pdb,irotam_pdb,
- & itorkcc,ithetkcc,iliptranpar
+ & itorkcc,ithetkcc,iliptranpar,itube
character*256 outname,intname,pdbname,mol2name,statname,intinname,
& entname,prefix,secpred,rest2name,qname,cartname,tmpdir,
& mremd_rst_name,curdir,pref_orig
character*256 bondname,thetname,rotname,torname,tordname,
& fouriername,elename,sidename,scpname,sccorname,patname,
& thetname_pdb,rotname_pdb,liptranname,thetkccname,
- & torkccname
+ & torkccname,tubename
common /parfiles/ bondname,thetname,rotname,torname,tordname,
& fouriername,elename,sidename,scpname,sccorname,patname,
& thetname_pdb,rotname_pdb,liptranname,thetkccname,
- & torkccname
+ & torkccname,tubename
character*3 pot
C-----------------------------------------------------------------------
C INP - main input file
weights_(17)=wbond
weights_(18)=scal14
weights_(21)=wsccor
+ weights_(22)=wtube
+
C FG Master broadcasts the WEIGHTS_ array
call MPI_Bcast(weights_(1),n_ene,
& MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
wbond=weights(17)
scal14=weights(18)
wsccor=weights(21)
+ wtube=weights(22)
endif
time_Bcast=time_Bcast+MPI_Wtime()-time00
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
endif
+ if (TUBElog.gt.0) then
+C print *,"just before call"
+ call calctube(Etube)
+ else
+ Etube=0.0d0
+ endif
+
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
energia(22)=eliptran
energia(23)=Eafmforce
energia(24)=ethetacnstr
+ energia(25)=Etube
c Here are the energies showed per procesor if the are more processors
c per molecule then we sum it up in sum_energy subroutine
c print *," Processor",myrank," calls SUM_ENERGY"
eliptran=energia(22)
Eafmforce=energia(23)
ethetacnstr=energia(24)
+ Etube=energia(25)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
- & +ethetacnstr
+ & +ethetacnstr+wtube*Etube
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
& +Eafmforce
- & +ethetacnstr
+ & +ethetacnstr+wtube*Etube
#endif
energia(0)=etot
c detecting NaNQ
& +wturn3*gshieldc_t3(j,i)
& +wturn4*gshieldc_t4(j,i)
& +wel_loc*gshieldc_ll(j,i)
+ & +wtube*gg_tube(j,i)
+
enddo
& +wcorr*gshieldc_ec(j,i)
& +wturn4*gshieldc_t4(j,i)
& +wel_loc*gshieldc_ll(j,i)
+ & +wtube*gg_tube(j,i)
+
enddo
& +wturn4*gshieldc_loc_t4(j,i)
& +wel_loc*gshieldc_ll(j,i)
& +wel_loc*gshieldc_loc_ll(j,i)
-
-
-
-
-
+ & +wtube*gg_tube(j,i)
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& +wturn4*gshieldc_loc_t4(j,i)
& +wel_loc*gshieldc_ll(j,i)
& +wel_loc*gshieldc_loc_ll(j,i)
-
-
-
+ & +wtube*gg_tube(j,i)
#endif
& +wturn3*gshieldx_t3(j,i)
& +wturn4*gshieldx_t4(j,i)
& +wel_loc*gshieldx_ll(j,i)
+ & +wtube*gg_tube_sc(j,i)
eliptran=energia(22)
Eafmforce=energia(23)
ethetacnstr=energia(24)
+ etube=energia(25)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,
& ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
+ & etube,wtube,
& etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/
+ & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
& ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
+ & etube,wtube,
& etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/
+ & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
& *fac_shield(i)*fac_shield(j)
eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
+C#ifdef NEWCORR
C Derivatives in theta
gloc(nphi+i,icg)=gloc(nphi+i,icg)
& +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
& +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
& *fac_shield(i)*fac_shield(j)
-
+C#endif
C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
C Derivatives in shield mode
enddo
return
end
+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 calctube(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)=(c(1,i)+c(1,i+1))/2.0d0-tubecenter(1)
+ vectube(2)=(c(2,i)+c(2,i+1))/2.0d0-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...
+C .or.(iti.eq.10)
+ & ) cycle
+ vectube(1)=c(1,i+nres)-tubecenter(1)
+ vectube(2)=c(2,i+nres)-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
+ 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
+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
+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
+ 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
gliptranx(j,i)=0.0d0
gradafm(j,i)=0.0d0
grad_shield(j,i)=0.0d0
+ gg_tube(j,i)=0.0d0
+ gg_tube_sc(j,i)=0.0d0
C grad_shield_side is Cbeta sidechain gradient
do kk=1,maxshieldlist
grad_shield_side(j,kk,i)=0.0d0
icsa_pdb=42
itorkcc=43
ithetkcc=44
+ itube=45
C Lipidic input file for parameters range 60-79
iliptranpar=60
C input file for transfer sidechain and peptide group inside the
C
read (itorp,*,end=113,err=113) ntortyp
read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
+ itortyp(ntyp1)=ntortyp
do iblock=1,2
do i=-ntyp,-1
itortyp(i)=-itortyp(-i)
c ee(1,2,i)=0.0d0
c ee(2,1,i)=ee(1,2,i)
enddo
-c lprint=.true.
+ lprint=.true.
if (lprint) then
do i=1,nloctyp
write (iout,*) 'Type',i
enddo
enddo
endif
-c lprint=.false.
+ lprint=.false.
C
C Read electrostatic-interaction parameters
endif
enddo
enddo
+ write(iout,*) "tube param"
+ read(itube,*) epspeptube,sigmapeptube
+ sigmapeptube=sigmapeptube**6
+ sigeps=dsign(1.0D0,epspeptube)
+ 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
+ do i=1,ntyp
+ read(itube,*) epssctube,sigmasctube
+ sigmasctube=sigmasctube**6
+ sigeps=dsign(1.0D0,epssctube)
+ 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)
+ enddo
+
#ifdef OLDSCP
C
C Define the SC-p interaction constants (hard-coded; old style)
AFMlog=(index(controlcard,'AFM'))
selfguide=(index(controlcard,'SELFGUIDE'))
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
call reada(weightcard,'TEMP0',temp0,300.0d0)
call reada(weightcard,'WSHIELD',wshield,1.0d0)
call reada(weightcard,'WLT',wliptran,0.0D0)
+ call reada(weightcard,'WTUBE',wtube,1.0D0)
if (index(weightcard,'SOFT').gt.0) ipot=6
C 12/1/95 Added weight for the multi-body term WCORR
call reada(weightcard,'WCORRH',wcorr,1.0D0)
open (irotam_pdb,file=rotname_pdb,status='old',action='read')
#endif
#endif
+ call getenv_loc('TUBEPAR',tubename)
+ open (itube,file=tubename,status='old',readonly)
#ifndef OLDSCP
C
C 8/9/01 In the newest version SCp interaction constants are read from a file
projects / unres.git / commitdiff