From: Adam Sieradzan Date: Tue, 29 Sep 2015 17:54:56 +0000 (+0200) Subject: Merge branch 'lipid' into AFM X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=25618f9f83673a7063414fe1e17415d138f58da8;hp=3444e732000b8e1ae809ef2fab7a828697a9eb78 Merge branch 'lipid' into AFM Conflicts: source/wham/src-M/DIMENSIONS.ZSCOPT source/wham/src-M/energy_p_new.F source/wham/src-M/initialize_p.F source/wham/src-M/parmread.F --- diff --git a/PARAM/scinter_GB_ext_lip.parm b/PARAM/scinter_GB_ext_lip.parm index f4b894d..c508b3e 100644 --- a/PARAM/scinter_GB_ext_lip.parm +++ b/PARAM/scinter_GB_ext_lip.parm @@ -134,112 +134,111 @@ .0572167103 -.0468608825 .0151048455 .0084963678 .0278930397 .0076922911 .1033536738 -.0098256036 .0611385674 .0448303346 .0861379100 .0861379100 -C From here are lipid - 5.6053537261 6.2200032154 6.2159898496 6.4386968963 - 6.2098101231 5.9596374798 5.4310222662 4.8790085257 - 5.2638423822 5.4961954124 4.2928124162 4.3582250540 - 4.2667755921 3.6813845958 3.5405429647 3.7026985912 - 4.7177596156 3.3811992227 3.7936085458 4.5157490585 - 6.2200032154 6.2159898496 5.2638423822 5.2638423822 - - 6.6296746680 6.6717260508 6.8168876945 6.8378960152 - 6.3512039986 6.1425192423 5.5200206406 5.4858936024 - 4.9382573262 4.2795129166 4.0534857018 4.2086428525 - 3.5388086279 3.4209730327 2.7995744193 4.8164268553 - 3.8790755071 3.5878262177 4.6457432555 6.6296746680 - 6.6717260508 5.4858936024 5.4858936024 - - 6.6424306340 6.9715947250 6.9241225787 6.5291325353 - 6.1394821460 5.4415971840 5.2914044780 4.7881880526 - 4.1302408718 3.9405275117 3.9022731880 3.6078928145 - 3.0809713206 3.0182595342 4.7935507272 3.8711928134 - 3.6826168780 4.4872661030 6.6717260508 6.6424306340 - 5.2914044780 5.2914044780 - - 7.3271707934 7.2970176793 7.0494948235 6.2565687662 - 5.8133715113 5.8845382769 5.3463553241 4.7146928716 - 4.4864793359 4.3092352779 3.6928361281 4.0100566871 - 3.6584608489 4.5901676001 4.4681844010 4.6334334282 - 4.9023945082 6.8168876945 6.9715947250 5.8845382769 - 5.8845382769 - - 7.0106342493 6.8341607204 6.2047876082 5.5802122247 - 5.6716968726 5.2011344932 4.1966522349 4.2383931955 - 3.8954073453 3.5679086118 3.1136348127 2.7622561241 - 4.5092053374 4.0238412606 3.9799853578 4.8087180377 - 6.8378960152 6.9241225787 5.6716968726 5.6716968726 - - 6.5353606708 5.8844482033 5.1100976261 5.5515360893 - 4.9587559238 4.2030790774 4.0351600709 3.6506289806 - 3.4640416936 3.1611022087 2.7976192835 4.0422832763 - 3.4415952740 3.8620694983 4.5934246009 6.3512039986 - 6.5291325353 5.5515360893 5.5515360893 - - 5.2828811325 4.8298466315 4.7575327777 4.7249779042 - 3.5656384285 3.5613366567 3.9273045564 3.7664779118 - 3.5861699761 3.5523191627 4.6508672463 4.2424308260 - 4.0987727783 4.5079669929 6.1425192423 6.1394821460 - 4.7575327777 4.7575327777 - - 4.2222645754 4.1627149417 3.9786372270 3.1623319201 - 2.9379696734 3.1365985188 2.8025851476 2.6576343487 - 2.8733395837 3.8785891655 3.4443068252 3.4405132673 - 4.0915305763 5.5200206406 5.4415971840 4.1627149417 - 4.1627149417 - - 4.3860654351 3.6520014167 3.0264937881 2.5680579493 - 2.5001857370 2.1393094782 1.4707077389 1.7373181916 - 3.1045342752 2.0152329038 1.8196209136 3.6395852546 - 5.4858936024 5.2914044780 4.3860654351 4.3860654351 - - 2.5016557935 2.7879606493 2.1424756816 1.0741544747 - 1.3236438850 -.2712356678 .7613312181 2.8627348285 - 1.5002127649 -.2392577636 3.5381128985 4.9382573262 - 4.7881880526 3.6520014167 3.6520014167 - - 2.2903847447 2.2466771550 1.5053607856 1.4031797188 - .9676209687 1.3075541842 2.8595118604 1.9833053538 - -.0085773647 2.9557553666 4.2795129166 4.1302408718 - 3.0264937881 3.0264937881 - - 1.2800038243 .7689806100 .9119606366 -.0744371793 - .3669965427 2.4288563200 1.4975482174 -.2134101897 - 2.9415021479 4.0534857018 3.9405275117 2.5680579493 - 2.5680579493 - - -.6792428859 .4532383239 -.7590387660 -.3617034846 - 1.6803275058 .6775209988 -.5354837468 2.6208591363 - 4.2086428525 3.9022731880 2.5001857370 2.5001857370 - - .4136122500 -.1637152873 .2212564728 2.1287126020 - .3896388275 -.2354089650 2.3255326077 3.5388086279 - 3.6078928145 2.1393094782 2.1393094782 - - -3.5193531074 -2.0157307929 1.8355926160 2.7489174124 - 2.0098391558 1.7975718667 3.4209730327 3.0809713206 - 1.4707077389 1.4707077389 - - -1.3532688232 2.5504275249 2.8202873790 1.6253609671 - 1.8625091247 2.7995744193 3.0182595342 1.7373181916 - 1.7373181916 - - 3.7292778697 2.2944436481 -.0303565255 3.1115776177 - 4.8164268553 4.7935507272 3.1045342752 3.1045342752 - - -.0827362961 -1.6043113182 2.4439837435 3.8790755071 - 3.8711928134 2.0152329038 2.0152329038 - - -4.3897783280 2.3664634533 3.5878262177 3.6826168780 - 1.8196209136 1.8196209136 - - 4.1927969260 4.6457432555 4.4872661030 3.6395852546 - 3.6395852546 - - 6.6296746680 6.6717260508 5.4858936024 5.4858936024 - - 6.6424306340 5.2914044780 5.2914044780 - - 4.3860654351 4.3860654351 - - 4.3860654351 - + 2.252000 2.758289 2.828747 2.573807 + 2.573807 2.314337 3.329045 3.123838 + 1.679138 1.245463 2.417517 2.138542 + 2.902927 2.675567 2.916240 2.690005 + 2.934561 3.360153 2.904943 2.252400 + 2.822849 3.512073 2.021838 2.021838 + + 3.378400 3.464698 3.152444 3.152444 + 2.834641 4.077473 3.826132 2.056637 + 1.525464 2.961017 2.619323 3.555555 + 3.277081 3.571861 3.294765 3.594301 + 4.115573 3.558025 2.758779 3.457474 + 4.301648 2.476383 2.476383 + + 3.553200 3.232970 3.232970 2.907049 + 4.181627 3.923866 2.109172 1.564431 + 3.036654 2.686231 3.646378 3.360790 + 3.663100 3.378926 3.686114 4.220702 + 3.648911 2.829249 3.545792 4.411529 + 2.539639 2.539639 + + 2.941600 2.941600 2.645052 3.804760 + 3.570230 1.919084 1.423437 2.762977 + 2.444136 3.317750 3.057901 3.332965 + 3.074402 3.353905 3.840313 3.320055 + 2.574264 3.226230 4.013942 2.310756 + 2.310756 + + 2.941600 2.645052 3.804760 3.570230 + 1.919084 1.423437 2.762977 2.444136 + 3.317750 3.057901 3.332965 3.074402 + 3.353905 3.840313 3.320055 2.574264 + 3.226230 4.013942 2.310756 2.310756 + + 2.378400 3.421196 3.210309 1.725618 + 1.279938 2.484436 2.197739 2.983282 + 2.749629 2.996964 2.764467 3.015792 + 3.453165 2.985354 2.314748 2.900988 + 3.609290 2.077805 2.077805 + + 4.921200 4.617850 2.482205 1.841120 + 3.573723 3.161325 4.291286 3.955188 + 4.310965 3.976531 4.338049 4.967185 + 4.294266 3.329637 4.172909 5.191762 + 2.988806 2.988806 + + 4.333200 2.329199 1.727631 3.353434 + 2.966456 4.026765 3.711385 4.045231 + 3.731412 4.070646 4.661001 4.029562 + 3.124393 3.915686 4.871735 2.804572 + 2.804572 + + 1.252000 0.9286429 1.802551 1.594541 + 2.164482 1.994958 2.174409 2.005723 + 2.188069 2.505400 2.165986 1.679436 + 2.104775 2.618674 1.507524 1.507524 + + 0.6888000 1.337002 1.182715 1.605456 + 1.479715 1.612819 1.487700 1.622951 + 1.858324 1.606571 1.245684 1.561169 + 1.942343 1.118173 1.118173 + + 2.595200 2.295720 3.116286 2.872216 + 3.130577 2.887715 3.150245 3.607117 + 3.118450 2.417947 3.030323 3.770203 + 2.170439 2.170439 + + 2.030800 2.756675 2.540769 2.769316 + 2.554480 2.786715 3.190865 2.758589 + 2.138922 2.680631 3.335131 1.919976 + 1.919976 + + 3.742000 3.448923 3.759161 3.467535 + 3.782778 4.331384 3.744599 2.903442 + 3.638776 4.527216 2.606238 2.606238 + + 3.178800 3.464740 3.195954 3.486507 + 3.992146 3.451319 2.676042 3.353784 + 4.172640 2.402115 2.402115 + + -61.60000 -35.26000 3.800125 48.12500 + 35.00000 2.916757 3.655463 4.547977 + 2.618191 2.618191 + + -23.62500 3.505321 39.90000 28.44000 + 2.690483 3.371882 4.195157 2.415078 + 2.415078 + + 3.824000 4.378585 3.785405 2.935082 + 3.678429 4.576550 2.634640 2.634640 + + -1.447000 -28.07000 3.360749 4.211903 + 5.240276 3.016734 3.016734 + + -76.82000 2.905459 3.641304 4.530360 + 2.608049 2.608049 + + 2.252800 2.823350 3.512697 2.022197 + 2.022197 + + 3.538400 4.402332 2.534345 2.534345 + + 5.477200 3.153127 3.153127 + + 1.815200 1.815200 + + 1.815200 + diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index dc58df8..2635a2f 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -1671,11 +1671,11 @@ c Removing the velocity of the center of mass & "Time step reduced to",d_time, & " because of too large initial acceleration." endif - if(me.eq.king.or..not.out1file)then - write(iout,*) "Potential energy and its components" - call enerprint(potEcomp) +C if(me.eq.king.or..not.out1file)then +C write(iout,*) "Potential energy and its components" +C call enerprint(potEcomp) c write(iout,*) (potEcomp(i),i=0,n_ene) - endif +C endif potE=potEcomp(0)-potEcomp(20) totE=EK+potE itime=0 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 5402971..ee55c93 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2780,6 +2780,17 @@ c write(iout,*) 'b1=',b1(1,i-2) c write (iout,*) 'theta=', theta(i-1) enddo #else + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itortyp(itype(i-2)) + else + iti=ntortyp+1 + endif +c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itortyp(itype(i-1)) + else + iti1=ntortyp+1 + endif b1(1,i-2)=b(3,iti) b1(2,i-2)=b(5,iti) b2(1,i-2)=b(2,iti) @@ -2934,6 +2945,7 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo +C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) c write (iout,*) 'mu ',mu(:,i-2),i-2 cd write (iout,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) @@ -4151,7 +4163,7 @@ C Contribution to the local-electrostatic energy coming from the i-j pair & +a33*muij(4) c write (iout,*) 'i',i,' j',j,itype(i),itype(j), c & ' eel_loc_ij',eel_loc_ij -c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4) +C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) C Calculate patrial derivative for theta angle #ifdef NEWCORR geel_loc_ij=a22*gmuij1(1) diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index c881425..ae4d710 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -1130,11 +1130,12 @@ C---------------------- GB or BP potential ----------------------------- C now we start reading lipid do i=1,ntyp read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp) - print *,"WARNING!!" - do j=1,ntyp - epslip(i,j)=epslip(i,j)+0.05d0 - enddo +C print *,"WARNING!!" +C do j=1,ntyp +C epslip(i,j)=epslip(i,j)+0.05d0 +C enddo enddo + write(iout,*) epslip(1,1),"OK?" C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index d25c4eb..ebb59da 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -247,7 +247,7 @@ c Cutoff range for interactions bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick C endif - if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0)) + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick diff --git a/source/wham/src-M/COMMON.ALLPARM b/source/wham/src-M/COMMON.ALLPARM index 88b480f..5b6284a 100644 --- a/source/wham/src-M/COMMON.ALLPARM +++ b/source/wham/src-M/COMMON.ALLPARM @@ -66,8 +66,10 @@ & app_all(2,2,max_parm),bpp_all(2,2,max_parm), & ael6_all(2,2,max_parm),ael3_all(2,2,max_parm), & aad_all(ntyp,2,max_parm),bad_all(ntyp,2,max_parm), - & aa_all(ntyp,ntyp,max_parm),bb_all(ntyp,ntyp,max_parm), + & aa_aq_all(ntyp,ntyp,max_parm),bb_aq_all(ntyp,ntyp,max_parm), + & aa_lip_all(ntyp,ntyp,max_parm),bb_lip_all(ntyp,ntyp,max_parm), & augm_all(ntyp,ntyp,max_parm),eps_all(ntyp,ntyp,max_parm), + & epslip_all(ntyp,ntyp,max_parm), & sigma_all(ntyp,ntyp,max_parm),r0_all(ntyp,ntyp,max_parm), & chi_all(ntyp,ntyp,max_parm),chip_all(ntyp,max_parm), & alp_all(ntyp,max_parm),ebr_all(max_parm),d0cm_all(max_parm), @@ -98,7 +100,8 @@ & v0_all,v1_all,v2_all,vlor1_all,vlor2_all,vlor3_all,v1c_all, & v1s_all,v2c_all,v2s_all,b1_all,b2_all,cc_all,dd_all,ee_all, & ctilde_all,dtilde_all,b1tilde_all,app_all,bpp_all,ael6_all, - & ael3_all,aad_all,bad_all,aa_all,bb_all,augm_all, + & ael3_all,aad_all,bad_all,aa_aq_all,bb_aq_all,augm_all, + & aa_lip_all,bb_lip_all,epslip_all, & eps_all,sigma_all,r0_all,chi_all,chip_all,alp_all,ebr_all, & d0cm_all,akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all, & v1sccor_all,v2sccor_all,nbondterm_all, diff --git a/source/wham/src-M/COMMON.CHAIN b/source/wham/src-M/COMMON.CHAIN index b147f08..24b8c56 100644 --- a/source/wham/src-M/COMMON.CHAIN +++ b/source/wham/src-M/COMMON.CHAIN @@ -13,6 +13,8 @@ &nend_sup,chain_length,tabperm(maxperm,maxsym),nperm, & nstart_seq,ishift_pdb double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss, - &sssgrad - common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad + &sssgrad, + & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick + common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad, + & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick diff --git a/source/wham/src-M/COMMON.IOUNITS b/source/wham/src-M/COMMON.IOUNITS index 23783bb..188d55e 100644 --- a/source/wham/src-M/COMMON.IOUNITS +++ b/source/wham/src-M/COMMON.IOUNITS @@ -10,11 +10,12 @@ C----------------------------------------------------------------------- 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 + & istat,ientin,ientout,isidep1,ibond,ihist,izsc,idistr, + & iliptranpar 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 + & idistr,iliptranpar character*256 outname,intname,pdbname,mol2name,statname,intinname, & entname,restartname,prefix,scratchdir,sidepname,pdbfile, & histname,zscname @@ -23,9 +24,11 @@ C General I/O units & files & sidepname,pdbfile,histname,zscname C Parameter files character*256 bondname,thetname,rotname,torname,tordname, - & fouriername,elename,sidename,scpname,sccorname,patname + & fouriername,elename,sidename,scpname,sccorname,patname, + & liptranname common /parfiles/ thetname,rotname,torname,tordname,bondname, - & fouriername,elename,sidename,scpname,sccorname,patname + & fouriername,elename,sidename,scpname,sccorname,patname, + & liptranname character*3 pot C----------------------------------------------------------------------- C INP - main input file diff --git a/source/wham/src-M/enecalc1.F b/source/wham/src-M/enecalc1.F index fee94cf..a5d25b3 100644 --- a/source/wham/src-M/enecalc1.F +++ b/source/wham/src-M/enecalc1.F @@ -228,7 +228,7 @@ c call intout endif C write (iout,*) "Czy tu dochodze" potE(iii+1,iparm)=energia(0) - do k=1,21 + do k=1,22 enetb(k,iii+1,iparm)=energia(k) enddo #ifdef DEBUG diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index a7b0798..cba6b5e 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -91,6 +91,11 @@ C C 21/5/07 Calculate local sicdechain correlation energy C call eback_sc_corr(esccor) + + if (wliptran.gt.0) then + call Eliptransfer(eliptran) + endif + C C 12/1/95 Multi-body terms C @@ -116,6 +121,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 #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) @@ -125,6 +131,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 #endif energia(0)=etot energia(1)=evdw @@ -159,6 +166,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(20)=edihcnstr energia(21)=evdw_t energia(24)=ethetacnstr + energia(22)=eliptran c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -197,10 +205,12 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #else do i=1,nct @@ -216,10 +226,12 @@ C & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #endif enddo @@ -276,6 +288,7 @@ C------------------------------------------------------------------------ edihcnstr=energia(20) estr=energia(18) ethetacnstr=energia(24) + eliptran=energia(22) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -284,7 +297,8 @@ C------------------------------------------------------------------------ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),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,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss, + & eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -308,6 +322,7 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',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, @@ -316,7 +331,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,etot + & edihcnstr,ethetacnstr,ebr*nss,eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -339,6 +354,7 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -407,8 +423,8 @@ C Change 12/1/95 to calculate four-body interactions c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e1+e2 ij=icant(itypi,itypj) c ROZNICA z cluster @@ -422,7 +438,7 @@ cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@ -580,8 +596,8 @@ C rij=1.0D0/r_inv_ij r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e_augm+e1+e2 ij=icant(itypi,itypj) eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm) @@ -594,7 +610,7 @@ cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@ -726,8 +742,8 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives fac=(rrij*sigsq)**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -737,15 +753,15 @@ C to its derivatives eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux & /dabs(eps(itypi,itypj)) eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij endif if (calc_grad) then if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),15(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@ -822,6 +838,28 @@ C returning the ith atom to box if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize + if ((zi.gt.bordlipbot) + &.and.(zi.lt.bordliptop)) then +C the energy transfer exist + if (zi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) @@ -886,6 +924,33 @@ C returning jth atom to box if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize + if ((zj.gt.bordlipbot) + &.and.(zj.lt.bordliptop)) then +C the energy transfer exist + if (zj.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 + bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 +C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj) C checking the distance dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj @@ -931,6 +996,7 @@ c write (iout,*) i,j,xj,yj,zj if (sss.le.0.0) cycle C Calculate angle-dependent terms of energy and contributions to their C derivatives. + call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) @@ -944,13 +1010,13 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt - if (bb(itypi,itypj).gt.0) then + if (bb.gt.0) then evdw=evdw+evdwij*sss else evdw_t=evdw_t+evdwij*sss @@ -964,8 +1030,8 @@ c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj, c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)), c & aux*e2/eps(itypi,itypj) c if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa #ifdef DEBUG write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, @@ -1097,15 +1163,15 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm evdwij=evdwij*eps2rt*eps3rt - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij+e_augm else evdw_t=evdw_t+evdwij+e_augm @@ -1816,6 +1882,8 @@ c print *,"itilde3 i iti iti1",i,iti,iti1 do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) enddo +C write (iout,*) 'mumu',i,b1(1,iti),Ub2(1,i-2) + C Vectors and matrices dependent on a single virtual-bond dihedral. call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) @@ -2437,8 +2505,9 @@ C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) 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,'(a6,2i5,0pf7.3)') +C & 'eelloc',i,j,eel_loc_ij +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) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma @@ -2776,6 +2845,16 @@ C Cartesian derivatives enddo endif else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds + & .or.((i+5).gt.nres) + & .or.((i-1).le.0) +C end of changes suggested by Ana + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1 + & .or. itype(i+5).eq.ntyp1 + & .or. itype(i).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1) goto 178 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Fourth-order contributions @@ -2915,6 +2994,7 @@ C Remaining derivatives of this turn contribution gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) enddo endif + 178 continue endif return end @@ -7774,6 +7854,125 @@ cd write (2,*) 'eel_turn6',ekont*eel_turn6 return end crc------------------------------------------------- +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + subroutine Eliptransfer(eliptran) + 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.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' +C this is done by Adasko +C print *,"wchodze" +C structure of box: +C water +C--bordliptop-- buffore starts +C--bufliptop--- here true lipid starts +C lipid +C--buflipbot--- lipid ends buffore starts +C--bordlipbot--buffore ends + eliptran=0.0 + do i=1,nres +C do i=1,1 + if (itype(i).eq.ntyp1) cycle + + positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize)) + if (positi.le.0) positi=positi+boxzsize +C print *,i +C first for peptide groups +c for each residue check if it is in lipid or lipid water border area + if ((positi.gt.bordlipbot) + &.and.(positi.lt.bordliptop)) then +C the energy transfer exist + if (positi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 +C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 +C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran +C print *, "doing sscalefor top part" +C print *,i,sslip,fracinbuf,ssgradlip + else + eliptran=eliptran+pepliptran +C print *,"I am in true lipid" + endif +C else +C eliptran=elpitran+0.0 ! I am in water + endif + enddo +C print *, "nic nie bylo w lipidzie?" +C now multiply all by the peptide group transfer factor +C eliptran=eliptran*pepliptran +C now the same for side chains +CV do i=1,1 + do i=1,nres + if (itype(i).eq.ntyp1) cycle + 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) +C print *,positi,bordlipbot,buflipbot + if ((positi.gt.bordlipbot) + & .and.(positi.lt.bordliptop)) then +C the energy transfer exist + if (positi.lt.buflipbot) then + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0- + &((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + eliptran=eliptran+liptranene(itype(i)) +C print *,"I am in true lipid" + endif + endif ! if in lipid or buffor +C else +C eliptran=elpitran+0.0 ! I am in water + enddo + return + end + + +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE MATVEC2(A1,V1,V2) implicit real*8 (a-h,o-z) include 'DIMENSIONS' @@ -7937,4 +8136,34 @@ C----------------------------------------------------------------------- return end C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + double precision function sscalelip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscale=1.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscalelip=1.0d0+r*r*(2*r-3.0d0) +C else +C sscale=0d0 +C endif + return + end +C----------------------------------------------------------------------- + double precision function sscagradlip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscagrad=0.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscagradlip=r*(6*r-6.0d0) +C else +C sscagrad=0.0d0 +C endif + return + end + +C----------------------------------------------------------------------- diff --git a/source/wham/src-M/include_unres/COMMON.DERIV b/source/wham/src-M/include_unres/COMMON.DERIV index 79f8630..80cbed6 100644 --- a/source/wham/src-M/include_unres/COMMON.DERIV +++ b/source/wham/src-M/include_unres/COMMON.DERIV @@ -1,5 +1,6 @@ double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gvdwpp, & gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gvdwx,gradcorr,gradxorr, + & gliptranc,gliptranx, & gradcorr5,gradcorr6,gel_loc,gcorr3_turn,gcorr4_turn,gcorr6_turn, & gel_loc_loc,gel_loc_turn3,gel_loc_turn4,gel_loc_turn6,gcorr_loc, & g_corr5_loc,g_corr6_loc,gradb,gradbx,gsccorc,gsccorx,gsccor_loc, @@ -10,6 +11,8 @@ & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres), & gvdwc(3,maxres),gelc(3,maxres),gvdwpp(3,maxres), & gradx_scp(3,maxres), + & gliptranc(3,-1:maxres), + & gliptranx(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), diff --git a/source/wham/src-M/include_unres/COMMON.FFIELD b/source/wham/src-M/include_unres/COMMON.FFIELD index 0c169f7..6c432a9 100644 --- a/source/wham/src-M/include_unres/COMMON.FFIELD +++ b/source/wham/src-M/include_unres/COMMON.FFIELD @@ -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 + & r0_corr,wliptran 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,weights(max_ene), + & wturn6,wvdwpp,wbond,wliptran,weights(max_ene), & scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp common /potentials/ potname(5) character*3 potname diff --git a/source/wham/src-M/include_unres/COMMON.INTERACT b/source/wham/src-M/include_unres/COMMON.INTERACT index d38c470..7d6b59f 100644 --- a/source/wham/src-M/include_unres/COMMON.INTERACT +++ b/source/wham/src-M/include_unres/COMMON.INTERACT @@ -1,8 +1,10 @@ - double precision aa,bb,augm,aad,bad,app,bpp,ael6,ael3 + double precision aa_aq,bb_aq,augm,aad,bad,app,bpp,ael6,ael3, + & aa_lip,bb_lip 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 - common /interact/aa(ntyp,ntyp),bb(ntyp,ntyp),augm(ntyp,ntyp), + common /interact/aa_aq(ntyp,ntyp),bb_aq(ntyp,ntyp), + & augm(ntyp,ntyp),aa_lip(ntyp,ntyp),bb_lip(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), & iend(maxres,maxint_gr),itype(maxres),itel(maxres),itypro, @@ -12,8 +14,9 @@ C 12/1/95 Array EPS included in the COMMON block. double precision eps,sigma,sigmaii,rs0,chi,chip,chip0,alp,signa0, & sigii,sigma0,rr0,r0,r0e,r0d,rpp,epp,elpp6,elpp3,eps_scp,rscp, - & eps_orig + & eps_orig,epslip common /body/eps(ntyp,ntyp),sigma(ntyp,ntyp),sigmaii(ntyp,ntyp), + &epslip(ntyp,ntyp), & 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), @@ -26,4 +29,8 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters. & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp), & distchainmax,nbondterm(ntyp) &,vbldpDUM +C 01/29/15 Lipidic parameters + double precision pepliptran,liptranene + common /lipid/ pepliptran,liptranene(ntyp) + diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F index c66de63..9b03d11 100644 --- a/source/wham/src-M/initialize_p.F +++ b/source/wham/src-M/initialize_p.F @@ -62,6 +62,8 @@ C ihist=30 iweight=31 izsc=32 +C Lipidic input file for parameters range 60-79 + iliptranpar=60 C C Set default weights of the energy terms. C @@ -88,8 +90,10 @@ C enddo do i=1,ntyp do j=1,ntyp - aa(i,j)=0.0D0 - bb(i,j)=0.0D0 + aa_lip(i,j)=0.0D0 + bb_lip(i,j)=0.0D0 + aa_aq(i,j)=0.0D0 + bb_aq(i,j)=0.0D0 augm(i,j)=0.0D0 sigma(i,j)=0.0D0 r0(i,j)=0.0D0 diff --git a/source/wham/src-M/make_ensemble1.F b/source/wham/src-M/make_ensemble1.F index 5d7b750..699995e 100644 --- a/source/wham/src-M/make_ensemble1.F +++ b/source/wham/src-M/make_ensemble1.F @@ -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, + & escloc,eliptran, & 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 @@ -162,6 +162,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft estr=enetb(18,i,iparm) esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) + eliptran=enetb(22,i,iparm) #ifdef SPLITELE etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 @@ -171,7 +172,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 + & +wbond*estr+wliptran*eliptran #else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) @@ -181,7 +182,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 + & +wbond*estr+wliptran*eliptran #endif #ifdef MPI Fdimless(i)= diff --git a/source/wham/src-M/openunits.F b/source/wham/src-M/openunits.F index b9f54b7..5200b3e 100644 --- a/source/wham/src-M/openunits.F +++ b/source/wham/src-M/openunits.F @@ -48,6 +48,8 @@ C Get parameter filenames and open the parameter files. open (isidep,file=sidename,status='old') call mygetenv('SIDEP',sidepname) open (isidep1,file=sidepname,status="old") + call mygetenv('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index 0d61fee..0cf2755 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -142,6 +142,7 @@ c wsccor=ww(19) whpb=ww(15) wstrain=ww(15) + wliptran=ww(22) endif call card_concat(controlcard,.false.) @@ -251,6 +252,11 @@ c enddo enddo endif + read(iliptranpar,*) pepliptran + do i=1,ntyp + read(iliptranpar,*) liptranene(i) + enddo + close(iliptranpar) #ifdef CRYST_THETA C C Read the parameters of the probability distribution/energy expression @@ -1101,6 +1107,13 @@ C---------------------- GB or BP potential ----------------------------- read (isidep,*)(sigii(i),i=1,ntyp) read (isidep,*)(chip(i),i=1,ntyp) read (isidep,*)(alp(i),i=1,ntyp) + do i=1,ntyp + read (isidep,*)(epslip(i,j),j=i,ntyp) +C print *,"WARNING!!" +C do j=1,ntyp +C epslip(i,j)=epslip(i,j)+0.05d0 +C enddo + enddo C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -1139,6 +1152,7 @@ C Calculate the "working" parameters of SC interactions. do i=2,ntyp do j=1,i-1 eps(i,j)=eps(j,i) + epslip(i,j)=epslip(j,i) enddo enddo do i=1,ntyp @@ -1156,6 +1170,7 @@ C Calculate the "working" parameters of SC interactions. do i=1,ntyp do j=i,ntyp epsij=eps(i,j) + epsijlip=epslip(i,j) if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then rrij=sigma(i,j) else @@ -1167,10 +1182,16 @@ C Calculate the "working" parameters of SC interactions. epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) - aa(i,j)=epsij*rrij*rrij - bb(i,j)=-sigeps*epsij*rrij - aa(j,i)=aa(i,j) - bb(j,i)=bb(i,j) + aa_aq(i,j)=epsij*rrij*rrij + bb_aq(i,j)=-sigeps*epsij*rrij + aa_aq(j,i)=aa_aq(i,j) + bb_aq(j,i)=bb_aq(i,j) + sigeps=dsign(1.0D0,epsijlip) + epsijlip=dabs(epsijlip) + aa_lip(i,j)=epsijlip*rrij*rrij + bb_lip(i,j)=-sigeps*epsijlip*rrij + aa_lip(j,i)=aa_lip(i,j) + bb_lip(j,i)=bb_lip(i,j) if (ipot.gt.2) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F index 557d27b..684f090 100644 --- a/source/wham/src-M/readrtns.F +++ b/source/wham/src-M/readrtns.F @@ -82,6 +82,23 @@ c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick +C endif + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) + & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + if ((lipbufthick*2.0d0).gt.lipthick) + &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" + endif + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot call readi(controlcard,'SYM',symetr,1) write (iout,*) "DISTCHAINMAX",distchainmax write (iout,*) "delta",delta diff --git a/source/wham/src-M/store_parm.F b/source/wham/src-M/store_parm.F index 8c44422..f44b6f4 100644 --- a/source/wham/src-M/store_parm.F +++ b/source/wham/src-M/store_parm.F @@ -40,6 +40,7 @@ c Store weights ww_all(16,iparm)=wvdwpp ww_all(17,iparm)=wbond ww_all(19,iparm)=wsccor + ww_all(22,iparm)=wliptran c Store bond parameters vbldp0_all(iparm)=vbldp0 akp_all(iparm)=akp @@ -222,13 +223,16 @@ c Store the parameters of electrostatic interactions c Store sidechain parameters do i=1,ntyp do j=1,ntyp - aa_all(j,i,iparm)=aa(j,i) - bb_all(j,i,iparm)=bb(j,i) + aa_aq_all(j,i,iparm)=aa_aq(j,i) + bb_aq_all(j,i,iparm)=bb_aq(j,i) + aa_lip_all(j,i,iparm)=aa_lip(j,i) + bb_lip_all(j,i,iparm)=bb_lip(j,i) r0_all(j,i,iparm)=r0(j,i) sigma_all(j,i,iparm)=sigma(j,i) chi_all(j,i,iparm)=chi(j,i) augm_all(j,i,iparm)=augm(j,i) eps_all(j,i,iparm)=eps(j,i) + epslip_all(j,i,iparm)=epslip(j,i) enddo enddo do i=1,ntyp @@ -311,6 +315,7 @@ c Restore weights wvdwpp=ww_all(16,iparm) wbond=ww_all(17,iparm) wsccor=ww_all(19,iparm) + wliptran=ww_all(22,iparm) c Restore bond parameters vbldp0=vbldp0_all(iparm) akp=akp_all(iparm) @@ -492,13 +497,16 @@ c Restore the parameters of electrostatic interactions c Restore sidechain parameters do i=1,ntyp do j=1,ntyp - aa(j,i)=aa_all(j,i,iparm) - bb(j,i)=bb_all(j,i,iparm) + aa_aq(j,i)=aa_aq_all(j,i,iparm) + bb_aq(j,i)=bb_aq_all(j,i,iparm) + aa_lip(j,i)=aa_lip_all(j,i,iparm) + bb_lip(j,i)=bb_lip_all(j,i,iparm) r0(j,i)=r0_all(j,i,iparm) sigma(j,i)=sigma_all(j,i,iparm) chi(j,i)=chi_all(j,i,iparm) augm(j,i)=augm_all(j,i,iparm) eps(j,i)=eps_all(j,i,iparm) + epslip(j,i)=epslip_all(j,i,iparm) enddo enddo do i=1,ntyp diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index 15d6716..45573d1 100644 --- a/source/wham/src-M/wham_calc1.F +++ b/source/wham/src-M/wham_calc1.F @@ -84,7 +84,8 @@ c parameter (MaxHdim=200000) & eplus,eminus,logfac,tanhT,tt 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 + & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, + & eliptran integer ind_point(maxpoint),upindE,indE character*16 plik @@ -312,6 +313,8 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft estr=enetb(18,i,iparm) esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) + eliptran=enetb(22,i,iparm) + #ifdef DEBUG write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6), & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc, @@ -327,7 +330,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 + & +wbond*estr+wliptran*eliptran #else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) @@ -337,7 +340,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 + & +wbond*estr+wliptran*eliptran #endif #ifdef DEBUG write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3), @@ -776,7 +779,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 + & +wbond*estr+wliptran*eliptran eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ @@ -799,7 +802,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 + & +wbond*estr+wliptran*eliptran eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+