added nanostructures energy to wham, no differs
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 14 Mar 2017 17:12:48 +0000 (18:12 +0100)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 14 Mar 2017 17:12:48 +0000 (18:12 +0100)
16 files changed:
PARAM/tube_carbonano2.parm [new file with mode: 0644]
source/wham/src-M/COMMON.CHAIN
source/wham/src-M/COMMON.CONTROL
source/wham/src-M/COMMON.IOUNITS
source/wham/src-M/DIMENSIONS.ZSCOPT
source/wham/src-M/energy_p_new.F
source/wham/src-M/include_unres/COMMON.DERIV
source/wham/src-M/include_unres/COMMON.FFIELD
source/wham/src-M/include_unres/COMMON.INTERACT
source/wham/src-M/include_unres/COMMON.LOCAL
source/wham/src-M/initialize_p.F
source/wham/src-M/make_ensemble1.F
source/wham/src-M/openunits.F
source/wham/src-M/parmread.F
source/wham/src-M/readrtns.F
source/wham/src-M/wham_calc1.F

diff --git a/PARAM/tube_carbonano2.parm b/PARAM/tube_carbonano2.parm
new file mode 100644 (file)
index 0000000..ea98cb1
--- /dev/null
@@ -0,0 +1,25 @@
+1.1970470 5.3667307 0 0 0 0 3.0000000
+1.5539975 5.6438808 0 0 0 0 3.0000000
+1.6679316 5.6689787 0 0 0 0 3.0000000
+1.6606077 5.9381499 0 0 0 0 3.0000000
+1.7428987 5.8625088 0 0 0 0 3.0000000
+1.7310307 5.9950466 0 0 0 0 3.0000000
+1.6322831 5.8318806 0 0 0 0 3.0000000
+1.5348705 5.4955850 0 0 0 0 3.0000000
+1.3603992 5.3937664 0 0 0 0 3.0000000
+1.3228511 5.4371481 0 0 0 0 3.0000000
+1.1970470 5.3667307 0 0 0 0 3.0000000
+1.0325602 5.5439558 0 0 0 0 3.0000000
+0.98513186 5.3780737 0 0 0 0 3.0000000
+0.97556829 5.3995867 0 0 0 0 3.0000000
+0.90197319 5.4184709 0 0 0 0 3.0000000
+0.77024281 5.4679136 0 0 0 0 3.0000000
+0.75456488 5.4686551 0 0 0 0 3.0000000
+1.1983876 5.3466215 0 0 0 0 3.0000000
+0.96779823 5.2968884 0 0 0 0 3.0000000
+0.92065424 5.3752089 0 0 0 0 3.0000000
+1.1218165 5.6721835 0 0 0 0 3.0000000
+1.6679316 5.7029562 0 0 0 0 3.0000000
+1.6606077 5.9355397 0 0 0 0 3.0000000
+1.3228511 5.4343948 0 0 0 0 3.0000000
+1.3228511 5.4343948 0 0 0 0 3.0000000
index 24b8c56..7a4c4c6 100644 (file)
      & nstart_seq,ishift_pdb
       double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,
      &sssgrad,
-     & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+     & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick,
+     & tubecenter,tubeR0,
+     & buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick
       common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+      common /tube/ tubecenter(3),tubeR0,
+     & buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick
 
index fb20a3e..8e72c46 100644 (file)
@@ -1,5 +1,6 @@
       integer iscode,indpdb,outpdb,outmol2,icomparfunc,pdbint,
-     & ensembles,constr_dist,symetr,shield_mode,tor_mode
+     & ensembles,constr_dist,symetr,shield_mode,tor_mode,
+     & tubelog,genconstr
       logical refstr,pdbref,punch_dist,print_rms,caonly,verbose,
      & merge_helices,bxfile,cxfile,histfile,entfile,zscfile,
      & rmsrgymap,with_dihed_constr,check_conf,histout,with_theta_constr
@@ -8,4 +9,6 @@
      & merge_helices,bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap,
      & ensembles,with_dihed_constr,constr_dist,check_conf,histout,
      & with_theta_constr,
-     &symetr,tor_mode,shield_mode
+     &symetr,tor_mode,shield_mode,
+     & tubelog,genconstr
+
index 188d55e..c9c2d8c 100644 (file)
@@ -11,11 +11,11 @@ C General I/O units & files
       integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam,
      &        itorp,itordp,ifourier,ielep,isidep,iscpp,isccor,icbase,
      &        istat,ientin,ientout,isidep1,ibond,ihist,izsc,idistr,
-     &        iliptranpar
+     &        iliptranpar,itube
       common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,
      &        irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,isccor,
      &        icbase,istat,ientin,ientout,isidep1,ibond,ihist,izsc,
-     &        idistr,iliptranpar
+     &        idistr,iliptranpar,itube
       character*256 outname,intname,pdbname,mol2name,statname,intinname,
      &        entname,restartname,prefix,scratchdir,sidepname,pdbfile,
      &        histname,zscname
@@ -25,10 +25,10 @@ C General I/O units & files
 C Parameter files
       character*256 bondname,thetname,rotname,torname,tordname,
      &       fouriername,elename,sidename,scpname,sccorname,patname,
-     &       liptranname
+     &       liptranname,tubename
       common /parfiles/ thetname,rotname,torname,tordname,bondname,
      &       fouriername,elename,sidename,scpname,sccorname,patname,
-     &       liptranname
+     &       liptranname,tubename
       character*3 pot
 C-----------------------------------------------------------------------
 C INP    - main input file
index 6519938..a99d1ed 100644 (file)
@@ -3,7 +3,7 @@
 c Maximum number of structures in the database, energy components, proteins,
 c and structural classes
 c#ifdef JUBL
-      parameter (maxstr=200000,max_ene=25,maxprot=7,maxclass=5000)
+      parameter (maxstr=200000,max_ene=26,maxprot=7,maxclass=5000)
       parameter (maxclass1=10)
 c Maximum number of structures to be dealt with by one processor
       parameter (maxstr_proc=10000)
index 6c943cc..81df835 100644 (file)
@@ -104,6 +104,18 @@ C
         call Eliptransfer(eliptran)
       endif
 
+      if (TUBElog.eq.1) then
+      print *,"just before call"
+        call calctube(Etube)
+       print *,"just after call",etube
+       elseif (TUBElog.eq.2) then
+        call calctube2(Etube)
+       elseif (TUBElog.eq.3) then
+        call calcnano(Etube)
+       else
+       Etube=0.0d0
+       endif
+
 C 
 C 12/1/95 Multi-body terms
 C
@@ -136,7 +148,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
      & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
-     & +wliptran*eliptran
+     & +wliptran*eliptran+wtube*Etube
       else
       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
      & +wvdwpp*evdw1
@@ -146,7 +158,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
      & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
-     & +wliptran*eliptran
+     & +wliptran*eliptran+wtube*Etube
       endif
 #else
       if (shield_mode.gt.0) then
@@ -158,7 +170,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
      & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
-     & +wliptran*eliptran
+     & +wliptran*eliptran+wtube*Etube
       else
       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
      & +welec*fact(1)*(ees+evdw1)
@@ -168,7 +180,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
      & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
-     & +wliptran*eliptran
+     & +wliptran*eliptran+wtube*Etube
       endif
 #endif
       energia(0)=etot
@@ -205,6 +217,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
       energia(21)=evdw_t
       energia(24)=ethetacnstr
       energia(22)=eliptran
+      energia(25)=Etube
 c detecting NaNQ
 #ifdef ISNAN
 #ifdef AIX
@@ -260,6 +273,8 @@ C
      &                 +wturn4*gshieldc_loc_t4(j,i)
      &                 +wel_loc*gshieldc_ll(j,i)
      &                 +wel_loc*gshieldc_loc_ll(j,i)
+     &                +wtube*gg_tube(j,i)
+
 
           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
@@ -297,6 +312,8 @@ C
      &                 +wturn4*gshieldc_loc_t4(j,i)
      &                 +wel_loc*gshieldc_ll(j,i)
      &                 +wel_loc*gshieldc_loc_ll(j,i)
+     &                +wtube*gg_tube(j,i)
+
 
           gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)
      &                 +fact(1)*wscp*gradx_scp(j,i)+
@@ -350,6 +367,8 @@ C
      &                 +wturn3*gshieldx_t3(j,i)
      &                 +wturn4*gshieldx_t4(j,i)
      &                 +wel_loc*gshieldx_ll(j,i)
+     &                 +wtube*gg_tube_sc(j,i)
+
 
               else
           gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+
@@ -387,6 +406,8 @@ C
      &                 +wturn3*gshieldx_t3(j,i)
      &                 +wturn4*gshieldx_t4(j,i)
      &                 +wel_loc*gshieldx_ll(j,i)
+     &                 +wtube*gg_tube_sc(j,i)
+
 
          endif
         enddo
@@ -446,6 +467,7 @@ C------------------------------------------------------------------------
       estr=energia(18)
       ethetacnstr=energia(24)
       eliptran=energia(22)
+      Etube=energia(25)
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
      &  wvdwpp,
@@ -455,7 +477,7 @@ C------------------------------------------------------------------------
      &  eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
      &  eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
      &  esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,
-     & eliptran,wliptran,etot
+     & eliptran,wliptran,etube,wtube ,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -480,6 +502,7 @@ C------------------------------------------------------------------------
      & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ 
      & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
+     & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
      & 'ETOT=  ',1pE16.6,' (total)')
 #else
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond,
@@ -488,7 +511,7 @@ C------------------------------------------------------------------------
      &  ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
      &  eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
      &  eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
-     &  edihcnstr,ethetacnstr,ebr*nss,eliptran,wliptran,etot
+     &  edihcnstr,ethetacnstr,ebr*nss,eliptran,wliptran,etube,wtube,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -9207,7 +9230,6 @@ C simple Kihara potential
       include 'COMMON.LOCAL'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
-      include 'COMMON.NAMES'
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
@@ -9380,7 +9402,6 @@ C simple Kihara potential
       include 'COMMON.LOCAL'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
-      include 'COMMON.NAMES'
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
@@ -9642,7 +9663,6 @@ C simple Kihara potential
       include 'COMMON.LOCAL'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
-      include 'COMMON.NAMES'
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
index 7f8ddfb..aad4544 100644 (file)
@@ -8,7 +8,7 @@
      & 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_ll, gshieldc_loc_ll
+     & gshieldc_ll, gshieldc_loc_ll,gg_tube,gg_tube_sc
 
       integer nfl,icg
       logical calc_grad
@@ -27,7 +27,8 @@
      & gshieldx_t4(3,-1:maxres), gshieldc_t4(3,-1:maxres),
      & gshieldc_loc_t4(3,-1:maxres),
      & gshieldx_ll(3,-1:maxres), gshieldc_ll(3,-1:maxres),
-     & gshieldc_loc_ll(3,-1:maxres),
+     & gshieldc_loc_ll(3,-1:maxres),gg_tube(3,-1:maxres),
+     & gg_tube_sc(3,-1:maxres),
      & gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres),
      & gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres),
      & gradcorr5(3,maxres),gradcorr6(3,maxres),
index 6c432a9..0ddbf5d 100644 (file)
@@ -6,11 +6,11 @@ C-----------------------------------------------------------------------
       double precision wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
      &    wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
      &    wturn6,wvdwpp,wbond,weights,scal14,cutoff_corr,delt_corr,
-     &    r0_corr,wliptran
+     &    r0_corr,wliptran,wtube
       integer ipot,n_ene_comp
       common /ffield/ wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
      &    wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
-     &    wturn6,wvdwpp,wbond,wliptran,weights(max_ene),
+     &    wturn6,wvdwpp,wbond,wliptran,weights(max_ene),wtube,
      &    scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp
       common /potentials/ potname(5)
       character*3 potname
index 7d6b59f..7f6e1cc 100644 (file)
@@ -1,5 +1,8 @@
       double precision aa_aq,bb_aq,augm,aad,bad,app,bpp,ael6,ael3,
-     & aa_lip,bb_lip
+     & aa_lip,bb_lip,sc_aa_tube_par,sc_bb_tube_par,
+     & pep_aa_tube,pep_bb_tube,dcavtub,acavtub,bcavtub,ccavtub,
+     & dcavtubpep,acavtubpep,bcavtubpep,ccavtubpep
+
       integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro,ielstart,
      & ielend,nscp_gr,iscpstart,iscpend,iatsc_s,iatsc_e,iatel_s,
      & iatel_e,iatscp_s,iatscp_e,ispp,iscp,expon,expon2
@@ -17,6 +20,10 @@ C 12/1/95 Array EPS included in the COMMON block.
      & eps_orig,epslip
       common /body/eps(ntyp,ntyp),sigma(ntyp,ntyp),sigmaii(ntyp,ntyp),
      &epslip(ntyp,ntyp),
+     & sc_aa_tube_par(ntyp),sc_bb_tube_par(ntyp),
+     & pep_aa_tube,pep_bb_tube,acavtub(ntyp),dcavtub(ntyp),
+     & bcavtub(ntyp),ccavtub(ntyp),acavtubpep,dcavtubpep,
+     & bcavtubpep,ccavtubpep,
      & rs0(ntyp,ntyp),chi(ntyp,ntyp),chip(ntyp),chip0(ntyp),alp(ntyp),
      & sigma0(ntyp),sigii(ntyp),rr0(ntyp),r0(ntyp,ntyp),r0e(ntyp,ntyp),
      & r0d(ntyp,2),rpp(2,2),epp(2,2),elpp6(2,2),elpp3(2,2),
@@ -30,7 +37,10 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters.
      & distchainmax,nbondterm(ntyp)
      &,vbldpDUM
 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
 
index 3e68e82..2b3be2a 100644 (file)
@@ -42,12 +42,14 @@ C Virtual-bond lenghts
      & iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
      & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
      & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
-     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end
+     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,itube_start,
+     & itube_end
       common /peptbond/ vbl,vblinv,vblinv2,vbl_cis,vbl0
       common /indices/ loc_start,loc_end,ithet_start,ithet_end,
      & iphi_start,iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
      & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
      & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
-     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end
+     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,itube_start,
+     & itube_end
 C Inverses of the actual virtual bond lengths
       common /invlen/ vbld_inv(maxres2)
index 5800580..b4598e1 100644 (file)
@@ -62,6 +62,7 @@ C
       ihist=30
       iweight=31
       izsc=32
+      itube=45
 C Lipidic input file for parameters range 60-79
       iliptranpar=60
 C
@@ -268,18 +269,18 @@ c-------------------------------------------------------------------------
      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
      &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
-     &   "EAFM","ETHETC","EMPTY"/
+     &   "EAFM","ETHETC","EMPTY","ETUBE"/
       data wname /
      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
      &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
-     &   "WLIPTRAN","WAFM","WTHETC","WSHIELD"/
+     &   "WLIPTRAN","WAFM","WTHETC","WSHIELD","WTUBE"/
       data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
      &    1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
-     &    0.0d0,0.0,0.0d0,0.0d0,0.0d0,0.0d0/
+     &    0.0d0,0.0,0.0d0,0.0d0,0.0d0,0.0d0,0.0d0/
       data nprint_ene /22/
       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
-     &  16,15,17,20,21,24,22,23,1/
+     &  16,15,17,20,21,24,22,23,1,25/
       end 
 c---------------------------------------------------------------------------
       subroutine init_int_table
@@ -527,6 +528,9 @@ C Partition local interactions
       call int_bounds(nres-3,itau_start,itau_end)
       itau_start=itau_start+3
       itau_end=itau_end+3
+      call int_bounds(nres-1,itube_start,itube_end)
+      itube_start=itube_start
+      itube_end=itube_end
       if (lprint) then 
         write (iout,*) 'Processor:',MyID,
      & ' loc_start',loc_start,' loc_end',loc_end,
@@ -552,6 +556,8 @@ C Partition local interactions
       iphi_end=nct
       itau_start=4
       itau_end=nres
+      itube_start=1
+      itube_end=nres-1
 #endif
       return
       end 
index 80538f6..403e490 100644 (file)
@@ -23,7 +23,7 @@
       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
-     &      escloc,eliptran,
+     &      escloc,eliptran,etube,
      &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
      &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
       integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
@@ -163,6 +163,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             esccor=enetb(19,i,iparm)
             edihcnstr=enetb(20,i,iparm)
             eliptran=enetb(22,i,iparm)
+            etube=enetb(25,i,iparm)
+
 #ifdef SPLITELE
       if (shield_mode.gt.0) then
             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
@@ -174,7 +176,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*etube
       else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
      &      +wvdwpp*evdw1
@@ -184,7 +186,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*etube
       endif
 #else
       if (shield_mode.gt.0) then
@@ -196,7 +198,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*etube
       else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
      &      +ft(1)*welec*(ees+evdw1)
@@ -206,7 +208,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*etube
       endif
 
 #endif
index 5200b3e..85830ec 100644 (file)
@@ -50,6 +50,9 @@ C Get parameter filenames and open the parameter files.
       open (isidep1,file=sidepname,status="old")
       call mygetenv('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
+      call mygetenv('TUBEPAR',tubename)
+      open (itube,file=tubename,status='old')
+
 #ifndef OLDSCP
 C
 C 8/9/01 In the newest version SCp interaction constants are read from a file
index cd2bce9..2fc37d3 100644 (file)
@@ -28,10 +28,10 @@ C
       character*1 toronelet(-2:2) /"p","a","G","A","P"/
       logical lprint
       dimension blower(3,3,maxlob)
-      character*800 controlcard
+      character*900 controlcard
       character*256 bondname_t,thetname_t,rotname_t,torname_t,
      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
-     & sccorname_t
+     & sccorname_t,tubename_t
       integer ilen
       external ilen
       character*16 key
@@ -146,6 +146,7 @@ c
       wstrain=ww(15)
       wliptran=ww(22)
       wshield=ww(25)
+      wtube=ww(26)
       endif
 
       call card_concat(controlcard,.false.)
@@ -184,6 +185,12 @@ c Return if not own parameters
       call reads(controlcard,"SCPPAR",scpname_t,scpname)
       open (iscpp,file=scpname_t,status='old')
       rewind(iscpp)
+      call reads(controlcard,"TUBEPAR",tubename_t,tubename)
+      write(iout,*) tubename_t
+      write(iout,*) tubename
+      open (itube,file=tubename_t,status='old')
+      rewind(itube)
+
       write (iout,*) "Parameter set:",iparm
       write (iout,*) "Energy-term weights:"
       do i=1,n_ene
@@ -1233,6 +1240,26 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
          endif
         enddo
       enddo
+      write(iout,*) "tube param"
+      read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
+     & ccavtubpep,dcavtubpep,tubetranenepep
+      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,tubetranenepep
+      do i=1,ntyp
+       read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
+     & ccavtub(i),dcavtub(i),tubetranene(i)
+       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),tubetranene(i)
+      enddo
+
 C
 C Define the SC-p interaction constants
 C
index 06c7af5..ff5b18e 100644 (file)
        call reada(controlcard,'BOXX',boxxsize,100.0d0)
        call reada(controlcard,'BOXY',boxysize,100.0d0)
        call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+      call readi(controlcard,'TUBEMOD',tubelog,0)
+      write (iout,*) TUBElog,"TUBEMODE"
+      call readi(controlcard,'GENCONSTR',genconstr,0)
+
 c Cutoff range for interactions
       call reada(controlcard,"R_CUT",r_cut,15.0d0)
       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
@@ -106,6 +110,17 @@ C      endif
       write (iout,*) "einicheck",einicheck
       write (iout,*) "rescale_mode",rescale_mode
       call flush(iout)
+      if (TUBElog.gt.0) then
+       call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+       call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
+       call reada(controlcard,"ZTUBE",tubecenter(3),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
       bxfile=index(controlcard,"BXFILE").gt.0
       cxfile=index(controlcard,"CXFILE").gt.0
       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
index a3d1658..7c8fe50 100644 (file)
@@ -89,7 +89,7 @@ c      parameter (MaxHdim=200000)
       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
-     &  eliptran
+     &  eliptran,etube
 
       integer ind_point(maxpoint),upindE,indE
       character*16 plik
@@ -325,6 +325,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             esccor=enetb(19,i,iparm)
             edihcnstr=enetb(20,i,iparm)
             eliptran=enetb(22,i,iparm)
+            etube=enetb(25,i,iparm)
+
 
 #ifdef DEBUG
             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
@@ -343,7 +345,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
              else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
      &      +wvdwpp*evdw1
@@ -353,7 +355,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
              endif
 #else
       if (shield_mode.gt.0) then
@@ -365,7 +367,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
      &      +ft(1)*welec*(ees+evdw1)
@@ -375,7 +377,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
            endif
 
 #endif
@@ -693,6 +695,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
           edihcnstr=enetb(20,i,iparm)
 C          edihcnstr=0.0d0
           eliptran=enetb(22,i,iparm)
+            etube=enetb(25,i,iparm)
+
 #ifdef SPLITELE
       if (shield_mode.gt.0) then
             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
@@ -704,7 +708,7 @@ C          edihcnstr=0.0d0
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
         else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
      &      +wvdwpp*evdw1
@@ -714,7 +718,7 @@ C          edihcnstr=0.0d0
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
       endif
 #else
       if (shield_mode.gt.0) then
@@ -726,7 +730,7 @@ C          edihcnstr=0.0d0
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
        else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
      &      +ft(1)*welec*(ees+evdw1)
@@ -736,7 +740,7 @@ C          edihcnstr=0.0d0
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
        endif
 
 #endif
@@ -964,6 +968,7 @@ c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
           edihcnstr=enetb(20,t,iparm)
 C          edihcnstr=0.0d0
           eliptran=enetb(22,i,iparm)
+            etube=enetb(25,i,iparm)
 
           do k=0,nGridT
             betaT=startGridT+k*delta_T
@@ -1086,7 +1091,7 @@ c            write (iout,*) "ftbis",ftbis
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
 C     &            +ftprim(6)*evdw_t
      &            +ftprim(1)*wscp*evdw2
@@ -1117,7 +1122,7 @@ C     &            +ftprim(6)*evdw_t
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
      &            +ftprim(1)*wtor*etors+
      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
@@ -1142,7 +1147,7 @@ C     &            +ftprim(6)*evdw_t
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
      &           +ftprim(1)*welec*(ees+evdw1)
      &           +ftprim(1)*wtor*etors+
@@ -1169,7 +1174,7 @@ C     &            +ftprim(6)*evdw_t
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
      &           +ftprim(1)*wtor*etors+
      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+