From c1827efa66e69c93c6e2b2e7420b06b430c3550a Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 30 Mar 2016 22:43:39 +0200 Subject: [PATCH] introduction of infinite cylinder potential - currently without PBC --- PARAM/tube.parm | 26 +++++ source/unres/src_MD-M/COMMON.CHAIN | 4 +- source/unres/src_MD-M/COMMON.CONTROL | 4 +- source/unres/src_MD-M/COMMON.DERIV | 4 +- source/unres/src_MD-M/COMMON.FFIELD | 1 + source/unres/src_MD-M/COMMON.INTERACT | 6 +- source/unres/src_MD-M/COMMON.IOUNITS | 8 +- source/unres/src_MD-M/energy_p_new_barrier.F | 161 ++++++++++++++++++++++++-- source/unres/src_MD-M/gradient_p.F | 2 + source/unres/src_MD-M/initialize_p.F | 1 + source/unres/src_MD-M/parmread.F | 23 +++- source/unres/src_MD-M/readrtns_CSA.F | 10 ++ 12 files changed, 228 insertions(+), 22 deletions(-) create mode 100644 PARAM/tube.parm diff --git a/PARAM/tube.parm b/PARAM/tube.parm new file mode 100644 index 0000000..d1783aa --- /dev/null +++ b/PARAM/tube.parm @@ -0,0 +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 + diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index 84ebf1b..8443ef4 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -18,10 +18,12 @@ 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 diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index e13259d..88d7571 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -1,6 +1,6 @@ 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, @@ -11,7 +11,7 @@ & 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. diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV index 97395db..a830225 100644 --- a/source/unres/src_MD-M/COMMON.DERIV +++ b/source/unres/src_MD-M/COMMON.DERIV @@ -3,6 +3,7 @@ & 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, @@ -24,7 +25,8 @@ & 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), diff --git a/source/unres/src_MD-M/COMMON.FFIELD b/source/unres/src_MD-M/COMMON.FFIELD index a63fe78..c4f56bb 100644 --- a/source/unres/src_MD-M/COMMON.FFIELD +++ b/source/unres/src_MD-M/COMMON.FFIELD @@ -5,6 +5,7 @@ C 12/1/95 wcorr added 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, diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT index 448e829..c43ef6a 100644 --- a/source/unres/src_MD-M/COMMON.INTERACT +++ b/source/unres/src_MD-M/COMMON.INTERACT @@ -1,5 +1,7 @@ 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, @@ -8,6 +10,8 @@ & 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), diff --git a/source/unres/src_MD-M/COMMON.IOUNITS b/source/unres/src_MD-M/COMMON.IOUNITS index 28a14a1..9954d90 100644 --- a/source/unres/src_MD-M/COMMON.IOUNITS +++ b/source/unres/src_MD-M/COMMON.IOUNITS @@ -12,12 +12,12 @@ C General I/O units & files & 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 @@ -42,12 +42,12 @@ C Parameter files 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 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 4a47935..e86bb6e 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -55,6 +55,8 @@ C FG slaves as WEIGHTS array. 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) @@ -81,6 +83,7 @@ C FG slaves receive the WEIGHTS array 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 @@ -304,6 +307,13 @@ C print *,"za lipidami" 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 @@ -348,6 +358,7 @@ C 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" @@ -442,6 +453,7 @@ cMS$ATTRIBUTES C :: proc_proc 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 @@ -449,7 +461,7 @@ cMS$ATTRIBUTES C :: proc_proc & +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 @@ -458,7 +470,7 @@ cMS$ATTRIBUTES C :: proc_proc & +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 @@ -568,6 +580,8 @@ c enddo & +wturn3*gshieldc_t3(j,i) & +wturn4*gshieldc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) + & +wtube*gg_tube(j,i) + enddo @@ -591,6 +605,8 @@ c enddo & +wcorr*gshieldc_ec(j,i) & +wturn4*gshieldc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) + & +wtube*gg_tube(j,i) + enddo @@ -747,11 +763,7 @@ C print *,gradafm(1,13),"AFM" & +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)+ @@ -784,9 +796,7 @@ C print *,gradafm(1,13),"AFM" & +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 @@ -801,6 +811,7 @@ C print *,gradafm(1,13),"AFM" & +wturn3*gshieldx_t3(j,i) & +wturn4*gshieldx_t4(j,i) & +wel_loc*gshieldx_ll(j,i) + & +wtube*gg_tube_sc(j,i) @@ -1101,6 +1112,7 @@ C------------------------------------------------------------------------ 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, @@ -1109,6 +1121,7 @@ C------------------------------------------------------------------------ & 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)'/ @@ -1136,6 +1149,7 @@ C------------------------------------------------------------------------ & '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 @@ -1146,6 +1160,7 @@ C------------------------------------------------------------------------ & 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)'/ @@ -1172,6 +1187,7 @@ C------------------------------------------------------------------------ & '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 @@ -4682,6 +4698,7 @@ C fac_shield(j)=0.6 & *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 @@ -4689,7 +4706,7 @@ C Derivatives in theta 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 @@ -11670,4 +11687,126 @@ C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) 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 diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F index 4cfa18b..1d8f3d1 100644 --- a/source/unres/src_MD-M/gradient_p.F +++ b/source/unres/src_MD-M/gradient_p.F @@ -405,6 +405,8 @@ C end of zero grad for shielding 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 diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index bb87d16..6eafc7d 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -121,6 +121,7 @@ crc for ifc error 118 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 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 2d2c23b..a38147f 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -630,6 +630,7 @@ C Read torsional parameters 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) @@ -1132,7 +1133,7 @@ c ee(2,1,i)=0.0d0 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 @@ -1154,7 +1155,7 @@ c lprint=.true. enddo enddo endif -c lprint=.false. + lprint=.false. C C Read electrostatic-interaction parameters @@ -1359,6 +1360,24 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) 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) diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index c691135..89ddeac 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -145,6 +145,13 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword 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 @@ -645,6 +652,7 @@ C Read weights of the subsequent energy terms. 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) @@ -2131,6 +2139,8 @@ C Get parameter filenames and open the parameter files. 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 -- 1.7.9.5