From 5d299c1a16ab51f8206b8ee3b17c7bcabe9321b7 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 6 Jul 2017 10:59:09 +0200 Subject: [PATCH] working energy for shield and lipid wrong gradient --- PARAM/Lip_tran_initial_ext.parm | 25 + PARAM/scinter_GB_ext_lip.parm | 245 ++++++ source/unres/control.F90 | 6 + source/unres/data/calc_data.f90 | 2 +- source/unres/data/comm_local.f90 | 2 +- source/unres/data/control_data.f90 | 2 +- source/unres/data/energy_data.f90 | 22 +- source/unres/data/geometry_data.f90 | 4 + source/unres/data/io_units.f90 | 4 +- source/unres/energy.f90 | 1599 ++++++++++++++++++++++++++++++----- source/unres/io.f90 | 3 + source/unres/io_config.f90 | 121 ++- source/unres/minim.f90 | 9 +- 13 files changed, 1814 insertions(+), 230 deletions(-) create mode 100644 PARAM/Lip_tran_initial_ext.parm create mode 100644 PARAM/scinter_GB_ext_lip.parm diff --git a/PARAM/Lip_tran_initial_ext.parm b/PARAM/Lip_tran_initial_ext.parm new file mode 100644 index 0000000..65acddb --- /dev/null +++ b/PARAM/Lip_tran_initial_ext.parm @@ -0,0 +1,25 @@ + 0.02 + -2.09251 + -1.67129 + -2.43220 + -2.44579 + -2.30991 + -1.65770 + -3.05723 + -1.30442 + -0.42122 + 0.00000 + -0.35328 + 0.05435 + 0.29893 + 0.81526 + 0.86961 + 1.04625 + -0.17664 + 1.37236 + 1.34518 + -0.97831 + -1.67129 + -2.43220 + -0.35328 + -0.35328 diff --git a/PARAM/scinter_GB_ext_lip.parm b/PARAM/scinter_GB_ext_lip.parm new file mode 100644 index 0000000..6b31b91 --- /dev/null +++ b/PARAM/scinter_GB_ext_lip.parm @@ -0,0 +1,245 @@ +4 6 + 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.6748060017 2.7338810145 2.9664647229 2.8819636737 + 3.0210738150 2.8414286152 2.4773438660 2.4611943788 + 2.4653201213 2.4925087371 2.5734767751 2.4564026744 + 2.4847825281 2.4889289233 2.5089517645 2.5083338383 + 2.4220622723 2.2714609770 2.4520703089 2.7026129788 + 2.7338810145 2.9664647229 2.4653201213 2.4653201213 + + 4.9272154761 5.1054284230 4.2073516165 4.8513972837 + 2.7848875293 3.5829861634 7.8660217576 7.4299209847 + 1.9625939832 .7987769569 4.0580899681 1.8889021032 + 3.1987197026 3.2673274538 2.6848131904 2.0043027404 + 6.2446341910 8.1959452095 13.4748295858 2.6632376837 + 5.1054284230 4.2073516165 1.9625939832 1.9625939832 + + .8699023011 1.0540660014 .9385909298 1.0263274101 + 1.0835277045 1.0543183886 .7888686996 .8989305833 + 1.0039962875 1.2427518128 .8932801724 .9173928990 + 1.6157695657 1.4315860373 2.0498317879 1.4199615546 + .9933677971 1.4319625600 27.4951763288 .7788025286 + 1.0540660014 .9385909298 1.0039962875 1.0039962875 + + .0103697556 .0611385674 .0448303346 .0392831782 + .0854166338 .0398896619 .0249496569 .0232410908 + .0861379100 -.0754794185 -.0266146021 -.0163429099 + .0572167103 -.0468608825 .0151048455 .0084963678 + .0278930397 .0076922911 .1033536738 -.0098256036 + .0611385674 .0448303346 .0861379100 .0861379100 + + 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 + + -7.2938 -4.0239 3.800125 5.9101 + 3.8743 2.916757 3.655463 4.547977 + 2.618191 2.618191 + + -3.0057 3.505321 6.0636 3.5303 + 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 + + -0.1779 -3.4493 3.360749 4.211903 + 5.240276 3.016734 3.016734 + + -9.4378 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/control.F90 b/source/unres/control.F90 index 8d22bf0..6a8befd 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -229,6 +229,7 @@ iweight=31 izsc=32 #endif + iliptranpar=60 #if defined(WHAM_RUN) || defined(CLUSTER) ! ! setting the mpi variables for WHAM @@ -795,6 +796,9 @@ ! & " ivec_start",ivec_start," ivec_end",ivec_end iset_start=loc_start+2 iset_end=loc_end+2 + call int_bounds(nres,ilip_start,ilip_end) + ilip_start=ilip_start + ilip_end=ilip_end if (ndih_constr.eq.0) then idihconstr_start=1 idihconstr_end=0 @@ -1320,6 +1324,8 @@ iset_end=nres+1 iint_start=2 iint_end=nres-1 + ilip_start=1 + ilip_end=nres #endif !el common /przechowalnia/ ! deallocate(iturn3_start_all) diff --git a/source/unres/data/calc_data.f90 b/source/unres/data/calc_data.f90 index 3f40fd0..ed6d106 100644 --- a/source/unres/data/calc_data.f90 +++ b/source/unres/data/calc_data.f90 @@ -9,6 +9,6 @@ sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,& eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,& sigder,dsci_inv,dscj_inv - real(kind=8),dimension(3) :: erij,gg + real(kind=8),dimension(3) :: erij,gg,gg_lipi,gg_lipj !----------------------------------------------------------------------------- end module calc_data diff --git a/source/unres/data/comm_local.f90 b/source/unres/data/comm_local.f90 index ad29715..89a42b4 100644 --- a/source/unres/data/comm_local.f90 +++ b/source/unres/data/comm_local.f90 @@ -3,7 +3,7 @@ integer :: num_conti,j1,j2 real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,& - dz_normi,xmedi,ymedi,zmedi + dz_normi,xmedi,ymedi,zmedi,sslipi,sslipj,ssgradlipi,ssgradlipj real(kind=8),dimension(2,2) :: a_temp real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 diff --git a/source/unres/data/control_data.f90 b/source/unres/data/control_data.f90 index 414ac0c..ddad045 100644 --- a/source/unres/data/control_data.f90 +++ b/source/unres/data/control_data.f90 @@ -29,7 +29,7 @@ ! commom.control ! common /cntrl/ integer :: modecalc,iscode,indpdb,indback,indphi,iranconf,& - icheckgrad,iprint,i2ndstr,mucadyn,constr_dist,symetr + icheckgrad,iprint,i2ndstr,mucadyn,constr_dist,symetr,shield_mode logical :: minim,refstr,pdbref,overlapsc,& energy_dec,sideadd,lsecondary,read_cart,unres_pdb,& vdisulf,searchsc,lmuca,dccart,extconf,out1file,& diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index 1b05b85..65bf013 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -64,7 +64,7 @@ integer :: rescale_mode real(kind=8) :: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,& wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,& - wturn6,wvdwpp + wturn6,wvdwpp,wliptran,wshield,lipscale #ifdef CLUSTER real(kind=8) :: scalscp #endif @@ -88,7 +88,7 @@ !----------------------------------------------------------------------------- ! common.interact ! common /interact/ - real(kind=8),dimension(:,:),allocatable :: aa,bb,augm !(ntyp,ntyp) + real(kind=8),dimension(:,:),allocatable :: aa_aq,bb_aq,augm,aa_lip,bb_lip !(ntyp,ntyp) real(kind=8),dimension(:,:),allocatable :: aad,bad !(ntyp,2) real(kind=8),dimension(2,2) :: app,bpp,ael6,ael3 integer :: expon,expon2, nnt,nct,itypro @@ -101,7 +101,7 @@ ! 12/1/95 Array EPS included in the COMMON block. ! common /body/ real(kind=8),dimension(:,:),allocatable :: sigma !(0:ntyp1,0:ntyp1) - real(kind=8),dimension(:,:),allocatable :: eps,sigmaii,& + real(kind=8),dimension(:,:),allocatable :: eps,epslip,sigmaii,& rs0,chi,r0,r0e !(ntyp,ntyp) r0e !!! nie używane real(kind=8),dimension(:),allocatable :: chip,alp,sigma0,& sigii,rr0 !(ntyp) @@ -149,7 +149,8 @@ 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,& + ilip_start,ilip_end integer,dimension(:),allocatable :: ibond_displ,ibond_count,& ithet_displ,ithet_count,iphi_displ,iphi_count,iphi1_displ,& iphi1_count,ivec_displ,ivec_count,iset_displ,iset_count,& @@ -250,6 +251,17 @@ ! Varibles for cutoff on electorstatic real(kind=8) sss_ele_cut,sss_ele_grad integer xshift,yshift,zshift -!----------------------------------------------------------------------------- +!2 Jul 2017 lipidc parameters ----------------------------------------------------- + real(kind=8),dimension(:), allocatable :: liptranene + real(kind=8) :: pepliptran + +! 4 Jul 2017 parameters for shieliding + real(kind=8),dimension(:), allocatable :: long_r_sidechain, & + short_r_sidechain + real(kind=8) :: VSolvSphere,VSolvSphere_div,buff_shield + + + + end module energy_data diff --git a/source/unres/data/geometry_data.f90 b/source/unres/data/geometry_data.f90 index 022d4ab..f5b843b 100644 --- a/source/unres/data/geometry_data.f90 +++ b/source/unres/data/geometry_data.f90 @@ -21,6 +21,10 @@ integer :: nsup,nstart_sup,nstart_seq,chain_length,iprzes,nperm integer :: nend_sup,ishift_pdb !wham real(kind=8) :: rmssing,anatemp !wham + real(kind=8) :: buftubebot, buftubetop,bordtubebot,bordtubetop, & + tubebufthick + real(kind=8) :: buflipbot, bufliptop,bordlipbot,bordliptop, & + lipbufthick,lipthick integer,dimension(:,:),allocatable :: tabperm !(maxperm,maxsym) ! common /from_zscore/ in module.compare !----------------------------------------------------------------------------- diff --git a/source/unres/data/io_units.f90 b/source/unres/data/io_units.f90 index e470892..1a88d8c 100644 --- a/source/unres/data/io_units.f90 +++ b/source/unres/data/io_units.f90 @@ -14,7 +14,7 @@ integer :: 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 + ithep_pdb,irotam_pdb,iliptranpar,itube #ifdef WHAM_RUN ! el wham iounits integer :: isidep1,ihist,iweight,izsc,idistr @@ -43,7 +43,7 @@ ! common /parfiles/ character(len=256) :: bondname,thetname,rotname,torname,tordname,& fouriername,elename,sidename,scpname,sccorname,patname,& - thetname_pdb,rotname_pdb + thetname_pdb,rotname_pdb,liptranname !----------------------------------------------------------------------- ! INP - main input file ! IOUT - list file diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 1a69db0..014ce41 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -29,6 +29,8 @@ !----------------------------------------------------------------------------- ! Maximum number of SC local term fitting function coefficiants integer,parameter :: maxsccoef=65 +! Maximum number of local shielding effectors + integer,parameter :: maxcontsshi=50 !----------------------------------------------------------------------------- ! commom.calc common/calc/ !----------------------------------------------------------------------------- @@ -36,10 +38,13 @@ ! common /contacts/ ! Change 12/1/95 - common block CONTACTS1 included. ! common /contacts1/ + integer,dimension(:),allocatable :: num_cont !(maxres) integer,dimension(:,:),allocatable :: jcont !(maxconts,maxres) - real(kind=8),dimension(:,:),allocatable :: facont !(maxconts,maxres) + real(kind=8),dimension(:,:),allocatable :: facont,ees0plist !(maxconts,maxres) real(kind=8),dimension(:,:,:),allocatable :: gacont !(3,maxconts,maxres) + integer,dimension(:),allocatable :: ishield_list + integer,dimension(:,:),allocatable :: shield_list ! ! 12/26/95 - H-bonding contacts ! common /contacts_hb/ @@ -112,7 +117,13 @@ real(kind=8),dimension(:,:),allocatable :: gvdwc,gelc,gelc_long,& gvdwpp,gvdwc_scpp,gradx_scp,gvdwc_scp,ghpbx,ghpbc,& gradcorr,gradcorr_long,gradcorr5_long,gradcorr6_long,& - gcorr6_turn_long,gradxorr,gradcorr5,gradcorr6 !(3,maxres) + gcorr6_turn_long,gradxorr,gradcorr5,gradcorr6,gliptran,gliptranc,& + gliptranx, & + gshieldx,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,& + grad_shield !(3,maxres) ! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2) real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,& gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres) @@ -124,8 +135,11 @@ ! real(kind=8),dimension(:,:,:),allocatable :: dtheta !(3,2,maxres) real(kind=8),dimension(:,:),allocatable :: gscloc,gsclocx !(3,maxres) ! real(kind=8),dimension(:,:,:),allocatable :: dphi,dalpha,domega !(3,3,maxres) + real(kind=8),dimension(:,:,:),allocatable :: grad_shield_side, & + grad_shield_loc ! (3,maxcontsshileding,maxnres) ! integer :: nfl,icg ! common /deriv_loc/ + real(kind=8), dimension(:),allocatable :: fac_shield real(kind=8),dimension(3,5,2) :: derx,derx_turn ! common /deriv_scloc/ real(kind=8),dimension(:,:),allocatable :: dXX_C1tab,dYY_C1tab,& @@ -208,7 +222,7 @@ integer :: n_corr,n_corr1,ierror real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc - real(kind=8) :: eello_turn3,eello_turn4,estr,ebe + real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6 #ifdef MPI @@ -293,7 +307,7 @@ #endif ! ! Compute the side-chain and electrostatic interaction energy -! + print *, "Before EVDW" ! goto (101,102,103,104,105,106) ipot select case(ipot) ! Lennard-Jones potential. @@ -336,7 +350,11 @@ ! 50 continue end select ! continue - +! print *,"after EGB" +! shielding effect + if (shield_mode.eq.2) then + call set_shield_fac2 + endif !mc !mc Sep-06: egb takes care of dynamic ss bonds too !mc @@ -349,9 +367,10 @@ #ifdef TIMING time_vec=time_vec+MPI_Wtime()-time01 #endif -! print *,"Processor",myrank," left VEC_AND_DERIV" +! print *,"Processor",myrank," left VEC_AND_DERIV" if (ipot.lt.6) then #ifdef SPLITELE + print *,"after ipot if", ipot if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or. & wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0 & .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 & @@ -362,8 +381,9 @@ .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 & .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then #endif +! print *,"just befor eelec call" call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) -! write (iout,*) "ELEC calc" +! write (iout,*) "ELEC calc" else ees=0.0d0 evdw1=0.0d0 @@ -394,13 +414,13 @@ ! write (iout,*) "Soft-sphere SCP potential" call escp_soft_sphere(evdw2,evdw2_14) endif -!elwrite(iout,*) "in etotal before ebond",ipot +! write(iout,*) "in etotal before ebond",ipot ! ! Calculate the bond-stretching energy ! call ebond(estr) -!elwrite(iout,*) "in etotal afer ebond",ipot +! write(iout,*) "in etotal afer ebond",ipot ! ! Calculate the disulfide-bridge and other energy and the contributions @@ -491,7 +511,15 @@ Uconst=0.0d0 Uconst_back=0.0d0 endif -!elwrite(iout,*) "after Econstr" + call flush(iout) +! write(iout,*) "after Econstr" + + if (wliptran.gt.0) then +! print *,"PRZED WYWOLANIEM" + call Eliptransfer(eliptran) + else + eliptran=0.0d0 + endif #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 @@ -534,6 +562,7 @@ energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor + energia(22)=eliptran ! Here are the energies showed per procesor if the are more processors ! per molecule then we sum it up in sum_energy subroutine ! print *," Processor",myrank," calls SUM_ENERGY" @@ -574,7 +603,8 @@ logical :: reduce real(kind=8) :: evdw,evdw2,evdw2_14,ees,evdw1,ecorr,ecorr5,ecorr6 real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc - real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot + real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot, & + eliptran integer :: i #ifdef MPI integer :: ierr @@ -634,20 +664,21 @@ estr=energia(17) Uconst=energia(20) esccor=energia(21) + eliptran=energia(22) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & - +wbond*estr+Uconst+wsccor*esccor + +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & - +wbond*estr+Uconst+wsccor*esccor + +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran #endif energia(0)=etot ! detecting NaNQ @@ -770,7 +801,7 @@ !el local variables real(kind=8) :: etot,evdw,evdw2,ees,evdw1,ecorr,ecorr5,ecorr6,eel_loc real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc - real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor + real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran etot=energia(0) evdw=energia(1) @@ -800,6 +831,8 @@ estr=energia(17) Uconst=energia(20) esccor=energia(21) + eliptran=energia(22) + #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -808,7 +841,7 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,& edihcnstr,ebr*nss,& - Uconst,etot + Uconst,eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -832,6 +865,7 @@ 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & + 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/& 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,& @@ -840,7 +874,7 @@ ecorr,wcorr,& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,& - ebr*nss,Uconst,etot + ebr*nss,Uconst,eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -863,6 +897,7 @@ 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & + 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -887,7 +922,7 @@ ! include 'COMMON.NAMES' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTACTS' - real(kind=8),dimension(3) :: gg + real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj integer :: num_conti !el local variables integer :: i,itypi,iint,j,itypi1,itypj,k @@ -929,8 +964,8 @@ ! 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=e1+e2 !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -1051,7 +1086,7 @@ ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' - real(kind=8),dimension(3) :: gg + real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj logical :: scheck !el local variables integer :: i,iint,j,itypi,itypi1,k,itypj @@ -1084,8 +1119,8 @@ 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=e_augm+e1+e2 !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -1220,16 +1255,16 @@ ! Calculate whole angle-dependent part of epsilon and contributions ! to its derivatives fac=(rrij*sigsq)**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) !d write (iout,'(2(a3,i3,2x),15(0pf7.3))') !d & restyp(itypi),i,restyp(itypj),j, !d & epsi,sigm,chi1,chi2,chip1,chip2, @@ -1281,8 +1316,11 @@ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi real(kind=8) :: evdw,sig0ij real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& - dist_temp, dist_init + dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, & + sslipi,sslipj,faclip integer :: ii + real(kind=8) :: fracinbuf + !cccc energy_dec=.false. ! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 @@ -1290,6 +1328,7 @@ ! if (icall.eq.0) lprn=.false. !el ind=0 do i=iatsc_s,iatsc_e + print *,"I am in EVDW",i itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1303,6 +1342,29 @@ zi=dmod(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 + print *, sslipi,ssgradlipi dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1360,6 +1422,36 @@ if (yj.lt.0) yj=yj+boxysize zj=dmod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize +! print *,"tu",xi,yi,zi,xj,yj,zj +! print *,"tu2",j,j+nres,c(1,j),c(1,j+nres) +! this fragment set correct epsilon for lipid phase + 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 +!------------------------------------------------ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj @@ -1430,8 +1522,9 @@ !--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + faclip=fac + e1=fac*fac*aa!(itypi,itypj) + e2=fac*bb!(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -1441,8 +1534,8 @@ evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*sss_ele_cut 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!(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, & @@ -1451,8 +1544,9 @@ evdwij endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & - 'evdw',i,j,evdwij !,"egb" + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2e10.2,e11.3)')& + 'evdw',i,j,evdwij,xi,xj,rij !,"egb" +!C print *,i,j,c(1,i),c(1,j),c(2,i),c(2,j),c(3,i),c(3,j) ! if (energy_dec) write (iout,*) & ! 'evdw',i,j,evdwij @@ -1472,6 +1566,14 @@ gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac +!C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt)& + *(eps3rt*eps3rt)*sss_ele_cut/2.0d0*(faclip*faclip*& + (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))& + +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi + ! print *,'before sc_grad', gg(1),gg(2),gg(3) ! Calculate angular part of the gradient. call sc_grad @@ -1582,8 +1684,8 @@ !--------------------------------------------------------------- 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -1592,8 +1694,9 @@ evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij+e_augm 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_aq(itypi,itypj)/& + bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j,& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),& @@ -1639,7 +1742,7 @@ ! include 'COMMON.NAMES' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTACTS' - real(kind=8),dimension(3) :: gg + real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj !d print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct !el local variables integer :: i,iint,j,itypi,itypi1,itypj,k @@ -2098,7 +2201,7 @@ real(kind=8) :: auxvec(2),auxmat(2,2) integer :: i,iti1,iti,k,l real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2 - +! print *,"in set matrices" ! ! Compute the virtual-bond-torsional-angle dependent quantities needed ! to calculate the el-loc multibody terms of various order. @@ -2109,6 +2212,7 @@ #else do i=3,nres+1 #endif +! print *,i,"i" if (i .lt. nres+1) then sin1=dsin(phi(i)) cos1=dcos(phi(i)) @@ -2187,6 +2291,7 @@ else iti1=ntortyp+1 endif +! print *,iti,i,"iti",iti1,itype(i-1),itype(i-2) !d write (iout,*) '*******i',i,' iti1',iti !d write (iout,*) 'b1',b1(:,iti) !d write (iout,*) 'b2',b2(:,iti) @@ -2588,10 +2693,11 @@ !el local variables integer :: i,k,j real(kind=8) :: ees,evdw1,eel_loc,eello_turn3,eello_turn4 - real(kind=8) :: fac,t_eelecij + real(kind=8) :: fac,t_eelecij,fracinbuf !d write(iout,*) 'In EELEC' + print *,"IN EELEC" !d do i=1,nloctyp !d write(iout,*) 'Type',i !d write(iout,*) 'B1',B1(:,i) @@ -2630,6 +2736,8 @@ ! write (iout,*) 'i',i,' fac',fac enddo endif + print *,wel_loc,"wel_loc",wcorr4,wcorr5,wcorr6,wturn3,wturn4, & + wturn6 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or. & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then @@ -2637,11 +2745,15 @@ #ifdef TIMING time01=MPI_Wtime() #endif +! print *, "before set matrices" call set_matrices +! print *, "after set matrices" + #ifdef TIMING time_mat=time_mat+MPI_Wtime()-time01 #endif endif + print *, "after set matrices" !d do i=1,nres-1 !d write (iout,*) 'i=',i !d do k=1,3 @@ -2678,7 +2790,7 @@ ! - + print *,"before iturn3 loop" do i=iturn3_start,iturn3_end if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle @@ -2698,7 +2810,30 @@ zmedi=dmod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 - call eelecij(i,i+2,ees,evdw1,eel_loc) + if ((zmedi.gt.bordlipbot) & + .and.(zmedi.lt.bordliptop)) then +!C the energy transfer exist + if (zmedi.lt.buflipbot) then +!C what fraction I am in + fracinbuf=1.0d0- & + ((zmedi-bordlipbot)/lipbufthick) +!C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif + print *,i,sslipi,ssgradlipi + call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) num_cont_hb(i)=num_conti enddo @@ -2721,6 +2856,29 @@ if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=dmod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize + if ((zmedi.gt.bordlipbot) & + .and.(zmedi.lt.bordliptop)) then +!C the energy transfer exist + if (zmedi.lt.buflipbot) then +!C what fraction I am in + fracinbuf=1.0d0- & + ((zmedi-bordlipbot)/lipbufthick) +!C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif + num_conti=num_cont_hb(i) call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & @@ -2747,6 +2905,28 @@ if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=dmod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize + if ((zmedi.gt.bordlipbot) & + .and.(zmedi.lt.bordliptop)) then +!C the energy transfer exist + if (zmedi.lt.buflipbot) then +!C what fraction I am in + fracinbuf=1.0d0- & + ((zmedi-bordlipbot)/lipbufthick) +!C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif ! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) @@ -2795,8 +2975,8 @@ !el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 real(kind=8),dimension(4) :: muij real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& - dist_temp, dist_init - integer xshift,yshift,zshift + dist_temp, dist_init,rlocshield,fracinbuf + integer xshift,yshift,zshift,ilist,iresshield !el integer :: num_conti,j1,j2 !el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,& !el dz_normi,xmedi,ymedi,zmedi @@ -2861,6 +3041,29 @@ 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 + isubchap=0 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj @@ -2922,11 +3125,30 @@ evdwij=ev1+ev2 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac - eesij=el1+el2 +! eesij=el1+el2 + if (shield_mode.gt.0) then +!C fac_shield(i)=0.4 +!C fac_shield(j)=0.6 + el1=el1*fac_shield(i)**2*fac_shield(j)**2 + el2=el2*fac_shield(i)**2*fac_shield(j)**2 + eesij=(el1+el2) + ees=ees+eesij*sss_ele_cut +!C FOR NOW SHIELD IS NOT USED WITH LIPSCALE +!C & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + else + fac_shield(i)=1.0 + fac_shield(j)=1.0 + eesij=(el1+el2) + ees=ees+eesij & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)*sss_ele_cut +!C print *,"TUCC",(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 + endif + ! 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) - ees=ees+eesij*sss_ele_cut - evdw1=evdw1+evdwij*sss_ele_cut +! ees=ees+eesij*sss_ele_cut + evdw1=evdw1+evdwij*sss_ele_cut & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) !d write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') !d & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, !d & 1.0D0/dsqrt(rrmij),evdwij,eesij, @@ -2943,8 +3165,10 @@ ! Calculate contributions to the Cartesian gradient. ! #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut - facel=-3*rrmij*(el1+eesij)*sss_ele_cut + facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facel=-3*rrmij*(el1+eesij)*sss_ele_cut & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) fac1=fac erij(1)=xj*rmij erij(2)=yj*rmij @@ -2956,6 +3180,54 @@ ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & + (shield_mode.gt.0)) then +!C print *,i,j + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)& + *2.0*sss_ele_cut + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0 & + *sss_ele_cut + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) & + *2.0*sss_ele_cut + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 & + *sss_ele_cut + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield + enddo + enddo + do k=1,3 + gshieldc(k,i)=gshieldc(k,i)+ & + grad_shield(k,i)*eesij/fac_shield(i)*2.0 & + *sss_ele_cut + + gshieldc(k,j)=gshieldc(k,j)+ & + grad_shield(k,j)*eesij/fac_shield(j)*2.0 & + *sss_ele_cut + + gshieldc(k,i-1)=gshieldc(k,i-1)+ & + grad_shield(k,i)*eesij/fac_shield(i)*2.0 & + *sss_ele_cut + + gshieldc(k,j-1)=gshieldc(k,j-1)+ & + grad_shield(k,j)*eesij/fac_shield(j)*2.0 & + *sss_ele_cut + + enddo + endif + + ! do k=1,3 ! ghalf=0.5D0*ggg(k) ! gelc(k,i)=gelc(k,i)+ghalf @@ -2966,6 +3238,15 @@ gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo + gelc_long(3,j)=gelc_long(3,j)+ & + ssgradlipj*eesij/2.0d0*lipscale**2& + *sss_ele_cut + + gelc_long(3,i)=gelc_long(3,i)+ & + ssgradlipi*eesij/2.0d0*lipscale**2& + *sss_ele_cut + + ! ! Loop over residues i+1 thru j-1. ! @@ -2974,9 +3255,13 @@ !grad gelc(l,k)=gelc(l,k)+ggg(l) !grad enddo !grad enddo - ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj - ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj - ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj + ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ! do k=1,3 ! ghalf=0.5D0*ggg(k) ! gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -2987,8 +3272,13 @@ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo -! -! Loop over residues i+1 thru j-1. + +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ & + sss_ele_cut*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ & + sss_ele_cut*ssgradlipi*evdwij/2.0d0*lipscale**2 +!! Loop over residues i+1 thru j-1. ! !grad do k=i+1,j-1 !grad do l=1,3 @@ -2996,7 +3286,9 @@ !grad enddo !grad enddo #else - facvdw=(ev1+evdwij)*sss_ele_cut + facvdw=(ev1+evdwij)*sss_ele_cut & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facel=(el1+eesij)*sss_ele_cut fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) @@ -3028,13 +3320,22 @@ !grad enddo !grad enddo ! 9/28/08 AL Gradient compotents will be summed only at the end - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + ggg(1)=facvdw*xj & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(2)=facvdw*yj & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(3)=facvdw*zj & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ & + sss_ele_cut*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ & + sss_ele_cut*ssgradlipi*evdwij/2.0d0*lipscale**2 + #endif ! ! Angular part @@ -3051,7 +3352,10 @@ !d print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), !d & (dcosg(k),k=1,3) do k=1,3 - ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss_ele_cut + ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss_ele_cut & + *fac_shield(i)**2*fac_shield(j)**2 & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + enddo ! do k=1,3 ! ghalf=0.5D0*ggg(k) @@ -3071,11 +3375,17 @@ gelc(k,i)=gelc(k,i) & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)& - *sss_ele_cut + *sss_ele_cut & + *fac_shield(i)**2*fac_shield(j)**2 & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + gelc(k,j)=gelc(k,j) & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)& - *sss_ele_cut + *sss_ele_cut & + *fac_shield(i)**2*fac_shield(j)**2 & + *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -3267,6 +3577,68 @@ ! 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) + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 + endif + eel_loc_ij=eel_loc_ij & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) +!C Now derivative over eel_loc + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & + (shield_mode.gt.0)) then +!C print *,i,j + + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij & + /fac_shield(i)& + *sss_ele_cut + gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i) & + *sss_ele_cut + + gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)& + +rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij & + /fac_shield(j) & + *sss_ele_cut + gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j) & + *sss_ele_cut + + gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)& + +rlocshield + + enddo + enddo + + do k=1,3 + gshieldc_ll(k,i)=gshieldc_ll(k,i)+ & + grad_shield(k,i)*eel_loc_ij/fac_shield(i) & + *sss_ele_cut + gshieldc_ll(k,j)=gshieldc_ll(k,j)+ & + grad_shield(k,j)*eel_loc_ij/fac_shield(j) & + *sss_ele_cut + gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+ & + grad_shield(k,i)*eel_loc_ij/fac_shield(i) & + *sss_ele_cut + gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+ & + grad_shield(k,j)*eel_loc_ij/fac_shield(j) & + *sss_ele_cut + + enddo + endif + + ! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij ! eel_loc_ij=0.0 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & @@ -3281,11 +3653,16 @@ gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) & - *sss_ele_cut + *sss_ele_cut & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) & - *sss_ele_cut + *sss_ele_cut & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) ! Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) ! do l=1,3 ! ggg(1)=(agg(1,1)*muij(1)+ & @@ -3316,6 +3693,14 @@ !grad gel_loc(l,i)=gel_loc(l,i)+ghalf !grad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ & + ssgradlipj*eel_loc_ij/2.0d0*lipscale/ & + ((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + + gel_loc_long(3,i)=gel_loc_long(3,i)+ & + ssgradlipi*eel_loc_ij/2.0d0*lipscale/ & + ((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + !grad do k=i+1,j2 !grad do l=1,3 !grad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -3414,6 +3799,12 @@ else ees0pij=0 endif + if (shield_mode.eq.0) then + fac_shield(i)=1.0d0 + fac_shield(j)=1.0d0 + else + ees0plist(num_conti,i)=j + endif ! ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2 if (ees0tmp.gt.0) then @@ -3423,10 +3814,12 @@ endif ! ees0mij=0.0D0 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) & - *sss_ele_cut + *sss_ele_cut & + *fac_shield(i)*fac_shield(j) ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) & - *sss_ele_cut + *sss_ele_cut & + *fac_shield(i)*fac_shield(j) ! Diagnostics. Comment out or remove after debugging! ! ees0p(num_conti,i)=0.5D0*fac3*ees0pij @@ -3505,28 +3898,28 @@ gacontp_hb1(k,num_conti,i)= & !ghalfp+ (ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) & - *sss_ele_cut + *sss_ele_cut*fac_shield(i)*fac_shield(j) gacontp_hb2(k,num_conti,i)= & !ghalfp+ (ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)& - *sss_ele_cut + *sss_ele_cut*fac_shield(i)*fac_shield(j) gacontp_hb3(k,num_conti,i)=gggp(k) & - *sss_ele_cut + *sss_ele_cut*fac_shield(i)*fac_shield(j) gacontm_hb1(k,num_conti,i)= & !ghalfm+ (ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) & - *sss_ele_cut + *sss_ele_cut*fac_shield(i)*fac_shield(j) gacontm_hb2(k,num_conti,i)= & !ghalfm+ (ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) & - *sss_ele_cut + *sss_ele_cut*fac_shield(i)*fac_shield(j) gacontm_hb3(k,num_conti,i)=gggm(k) & - *sss_ele_cut + *sss_ele_cut*fac_shield(i)*fac_shield(j) enddo ! Diagnostics. Comment out or remove after debugging! @@ -3596,11 +3989,38 @@ !el dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,& !el num_conti,j1,j2 !el local variables - integer :: i,j,l - real(kind=8) :: eello_turn3 + integer :: i,j,l,k,ilist,iresshield + real(kind=8) :: eello_turn3,zj,fracinbuf,eello_t3, rlocshield j=i+2 ! write (iout,*) "eturn3",i,j,j1,j2 + zj=(c(3,j)+c(3,j+1))/2.0d0 + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + if ((zj.lt.0)) write (*,*) "CHUJ" + 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 + a_temp(1,1)=a22 a_temp(1,2)=a23 a_temp(2,1)=a32 @@ -3619,9 +4039,60 @@ call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1)) call transpose2(auxmat(1,1),auxmat1(1,1)) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) - eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) + if (shield_mode.eq.0) then + fac_shield(i)=1.0d0 + fac_shield(j)=1.0d0 + endif + + eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t3= & + 0.5d0*(pizda(1,1)+pizda(2,2)) & + *fac_shield(i)*fac_shield(j) + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & + (shield_mode.gt.0)) then +!C print *,i,j + + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i) + gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1) & + +rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j) + gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j) + gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1) & + +rlocshield + + enddo + enddo + + do k=1,3 + gshieldc_t3(k,i)=gshieldc_t3(k,i)+ & + grad_shield(k,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,j)=gshieldc_t3(k,j)+ & + grad_shield(k,j)*eello_t3/fac_shield(j) + gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+ & + grad_shield(k,i)*eello_t3/fac_shield(i) + gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+ & + grad_shield(k,j)*eello_t3/fac_shield(j) + enddo + endif + !d write (2,*) 'i,',i,' j',j,'eello_turn3', !d & 0.5d0*(pizda(1,1)+pizda(2,2)), !d & ' eello_turn3_num',4*eello_turn3_num @@ -3635,7 +4106,10 @@ call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & - +0.5d0*(pizda(1,1)+pizda(2,2)) + +0.5d0*(pizda(1,1)+pizda(2,2)) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + ! Cartesian derivatives do l=1,3 ! ghalf1=0.5d0*agg(l,1) @@ -3648,29 +4122,49 @@ a_temp(2,2)=aggi(l,4)!+ghalf4 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & - +0.5d0*(pizda(1,1)+pizda(2,2)) + +0.5d0*(pizda(1,1)+pizda(2,2)) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) a_temp(2,1)=aggi1(l,3)!+agg(l,3) a_temp(2,2)=aggi1(l,4)!+agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & - +0.5d0*(pizda(1,1)+pizda(2,2)) + +0.5d0*(pizda(1,1)+pizda(2,2)) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 a_temp(2,2)=aggj(l,4)!+ghalf4 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & - +0.5d0*(pizda(1,1)+pizda(2,2)) + +0.5d0*(pizda(1,1)+pizda(2,2)) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) a_temp(2,2)=aggj1(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & - +0.5d0*(pizda(1,1)+pizda(2,2)) + +0.5d0*(pizda(1,1)+pizda(2,2)) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo + gshieldc_t3(3,i)=gshieldc_t3(3,i)+ & + ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j)=gshieldc_t3(3,j)+ & + ssgradlipj*eello_t3/4.0d0*lipscale + gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+ & + ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+ & + ssgradlipj*eello_t3/4.0d0*lipscale + return end subroutine eturn3 !----------------------------------------------------------------------------- @@ -3705,8 +4199,9 @@ !el dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,& !el num_conti,j1,j2 !el local variables - integer :: i,j,iti1,iti2,iti3,l - real(kind=8) :: eello_turn4,s1,s2,s3 + integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield + real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,& + rlocshield j=i+3 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC @@ -3722,6 +4217,32 @@ !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !d call checkint_turn4(i,a_temp,eello_turn4_num) ! write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2 + zj=(c(3,j)+c(3,j+1))/2.0d0 + 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 + a_temp(1,1)=a22 a_temp(1,2)=a23 a_temp(2,1)=a32 @@ -3742,7 +4263,57 @@ call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) - eello_turn4=eello_turn4-(s1+s2+s3) + if (shield_mode.eq.0) then + fac_shield(i)=1.0 + fac_shield(j)=1.0 + endif + + eello_turn4=eello_turn4-(s1+s2+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t4=-(s1+s2+s3) & + *fac_shield(i)*fac_shield(j) +!C Now derivative over shield: + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & + (shield_mode.gt.0)) then +!C print *,i,j + + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i) + gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) & + +rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j) + gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ & + rlocshield & + +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j) + gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) & + +rlocshield + + enddo + enddo + + do k=1,3 + gshieldc_t4(k,i)=gshieldc_t4(k,i)+ & + grad_shield(k,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,j)=gshieldc_t4(k,j)+ & + grad_shield(k,j)*eello_t4/fac_shield(j) + gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+ & + grad_shield(k,i)*eello_t4/fac_shield(i) + gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+ & + grad_shield(k,j)*eello_t4/fac_shield(j) + enddo + endif + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn4',i,j,-(s1+s2+s3) !d write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), @@ -3754,7 +4325,10 @@ s1=scalar2(b1(1,iti2),auxvec(1)) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) - gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) + gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + ! Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -3762,7 +4336,10 @@ call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,1)) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) - gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) + gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + ! Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -3773,7 +4350,10 @@ call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1)) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) - gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) + gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + ! Cartesian derivatives ! Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -3810,7 +4390,11 @@ call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) - gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) + gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + + a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -3824,7 +4408,11 @@ call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) - gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) + gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + + a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -3838,7 +4426,11 @@ call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) - gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) + gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -3853,8 +4445,20 @@ call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ! write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 - gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) + gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) & + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + enddo + gshieldc_t4(3,i)=gshieldc_t4(3,i)+ & + ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j)=gshieldc_t4(3,j)+ & + ssgradlipj*eello_t4/4.0d0*lipscale + gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+ & + ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+ & + ssgradlipj*eello_t4/4.0d0*lipscale + return end subroutine eturn4 !----------------------------------------------------------------------------- @@ -6985,10 +7589,11 @@ real(kind=8),dimension(3) :: gx,gx1 logical :: lprn !el local variables - integer :: i,j,k,l,jj,kk,ll + integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,& ees0mkl,ees,coeffpees0pij,coeffmees0mij,& - coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl + coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, & + rlocshield lprn=.false. eij=facont_hb(jj,i) @@ -7063,6 +7668,80 @@ !grad enddo ! write (iout,*) "ehbcorr",ekont*ees ehbcorr=ekont*ees + if (shield_mode.gt.0) then + j=ees0plist(jj,i) + l=ees0plist(kk,k) +!C print *,i,j,fac_shield(i),fac_shield(j), +!C &fac_shield(k),fac_shield(l) + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & + (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i) + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ & + rlocshield & + +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) & + +rlocshield + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j) + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ & + rlocshield & + +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) & + +rlocshield + enddo + enddo + + do ilist=1,ishield_list(k) + iresshield=shield_list(ilist,k) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k) + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ & + rlocshield & + +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) & + +rlocshield + enddo + enddo + do ilist=1,ishield_list(l) + iresshield=shield_list(ilist,l) + do m=1,3 + rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l) + gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ & + rlocshield & + +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l) + gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) & + +rlocshield + enddo + enddo + do m=1,3 + gshieldc_ec(m,i)=gshieldc_ec(m,i)+ & + grad_shield(m,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,j)=gshieldc_ec(m,j)+ & + grad_shield(m,j)*ehbcorr/fac_shield(j) + gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+ & + grad_shield(m,i)*ehbcorr/fac_shield(i) + gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+ & + grad_shield(m,j)*ehbcorr/fac_shield(j) + + gshieldc_ec(m,k)=gshieldc_ec(m,k)+ & + grad_shield(m,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,l)=gshieldc_ec(m,l)+ & + grad_shield(m,l)*ehbcorr/fac_shield(l) + gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+ & + grad_shield(m,k)*ehbcorr/fac_shield(k) + gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+ & + grad_shield(m,l)*ehbcorr/fac_shield(l) + + enddo + endif + endif return end function ehbcorr #ifdef MOMENT @@ -9333,7 +10012,7 @@ #ifdef MPI include 'mpif.h' #endif - real(kind=8),dimension(3,nres) :: gradbufc,gradbufx,gradbufc_sum,& + real(kind=8),dimension(3,-1:nres) :: gradbufc,gradbufx,gradbufc_sum,& gloc_scbuf !(3,maxres) real(kind=8),dimension(4*nres) :: glocbuf !(4*maxres) @@ -9411,7 +10090,7 @@ call flush(iout) #endif #ifdef SPLITELE - do i=1,nct + do i=0,nct do j=1,3 gradbufc(j,i)=wsc*gvdwc(j,i)+ & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ & @@ -9421,11 +10100,19 @@ wcorr5*gradcorr5_long(j,i)+ & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & - wstrain*ghpbc(j,i) + wstrain*ghpbc(j,i) & + +wliptran*gliptranc(j,i) & + +welec*gshieldc(j,i) & + +wcorr*gshieldc_ec(j,i) & + +wturn3*gshieldc_t3(j,i)& + +wturn4*gshieldc_t4(j,i)& + +wel_loc*gshieldc_ll(j,i) + + enddo enddo #else - do i=1,nct + do i=0,nct do j=1,3 gradbufc(j,i)=wsc*gvdwc(j,i)+ & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ & @@ -9436,7 +10123,14 @@ wcorr5*gradcorr5_long(j,i)+ & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & - wstrain*ghpbc(j,i) + wstrain*ghpbc(j,i) & + +wliptran*gliptranc(j,i) & + +welec*gshieldc(j,i)& + +wcorr*gshieldc_ec(j,i) & + +wturn4*gshieldc_t4(j,i) & + +wel_loc*gshieldc_ll(j,i) + + enddo enddo #endif @@ -9450,7 +10144,7 @@ enddo call flush(iout) #endif - do i=1,nres + do i=0,nres do j=1,3 gradbufc_sum(j,i)=gradbufc(j,i) enddo @@ -9468,7 +10162,7 @@ #ifdef TIMING ! time_allreduce=time_allreduce+MPI_Wtime()-time00 #endif - do i=nnt,nres + do i=0,nres do k=1,3 gradbufc(k,i)=0.0d0 enddo @@ -9493,7 +10187,7 @@ do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,nnt,-1 + do i=nres-2,-1,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -9516,7 +10210,7 @@ call flush(iout) #endif !el#undef DEBUG - do i=1,nres + do i=-1,nres do j=1,3 gradbufc_sum(j,i)=gradbufc(j,i) gradbufc(j,i)=0.0d0 @@ -9525,7 +10219,7 @@ do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,nnt,-1 + do i=nres-2,-1,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -9559,7 +10253,7 @@ !el if (.not.allocated(gradx)) allocate(gradx(3,nres,2)) !(3,maxres,2) !el if (.not.allocated(gradc)) allocate(gradc(3,nres,2)) !(3,maxres,2) !el----------------- - do i=1,nct + do i=-1,nct do j=1,3 #ifdef SPLITELE gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & @@ -9579,7 +10273,19 @@ wcorr6*gradcorr6(j,i)+ & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & - +wscloc*gscloc(j,i) + +wscloc*gscloc(j,i) & + +wliptran*gliptranc(j,i) & + +welec*gshieldc(j,i) & + +welec*gshieldc_loc(j,i) & + +wcorr*gshieldc_ec(j,i) & + +wcorr*gshieldc_loc_ec(j,i) & + +wturn3*gshieldc_t3(j,i) & + +wturn3*gshieldc_loc_t3(j,i) & + +wturn4*gshieldc_t4(j,i) & + +wturn4*gshieldc_loc_t4(j,i) & + +wel_loc*gshieldc_ll(j,i) & + +wel_loc*gshieldc_loc_ll(j,i) + #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & @@ -9598,13 +10304,33 @@ wcorr6*gradcorr6(j,i)+ & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & - +wscloc*gscloc(j,i) + +wscloc*gscloc(j,i) & + +wliptran*gliptranc(j,i) & + +welec*gshieldc(j,i) & + +welec*gshieldc_loc(j,) & + +wcorr*gshieldc_ec(j,i) & + +wcorr*gshieldc_loc_ec(j,i) & + +wturn3*gshieldc_t3(j,i) & + +wturn3*gshieldc_loc_t3(j,i) & + +wturn4*gshieldc_t4(j,i) & + +wturn4*gshieldc_loc_t4(j,i) & + +wel_loc*gshieldc_ll(j,i) & + +wel_loc*gshieldc_loc_ll(j,i) + + #endif 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*gsccorx(j,i) & - +wscloc*gsclocx(j,i) + +wscloc*gsclocx(j,i) & + +wliptran*gliptranx(j,i) & + +welec*gshieldx(j,i) & + +wcorr*gshieldx_ec(j,i) & + +wturn3*gshieldx_t3(j,i) & + +wturn4*gshieldx_t4(j,i) & + +wel_loc*gshieldx_ll(j,i) + enddo enddo #ifdef DEBUG @@ -9837,12 +10563,12 @@ enddo ! write (iout,*) "gg",(gg(k),k=1,3) do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) & + gvdwx(k,i)=gvdwx(k,i)-gg(k) +gg_lipi(k)& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv & *sss_ele_cut - gvdwx(k,j)=gvdwx(k,j)+gg(k) & + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv & *sss_ele_cut @@ -9861,8 +10587,8 @@ !grad enddo !grad enddo do l=1,3 - gvdwc(l,i)=gvdwc(l,i)-gg(l) - gvdwc(l,j)=gvdwc(l,j)+gg(l) + gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) enddo return end subroutine sc_grad @@ -11234,6 +11960,18 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' endif return end function sscagrad_ele + real(kind=8) function sscalelip(r) + real(kind=8) r,gamm + sscalelip=1.0d0+r*r*(2.0d0*r-3.0d0) + return + end function sscalelip +!C----------------------------------------------------------------------- + real(kind=8) function sscagradlip(r) + real(kind=8) r,gamm + sscagradlip=r*(6.0d0*r-6.0d0) + return + end function sscagradlip + !!!!!!!!!!!!!!! !----------------------------------------------------------------------------- subroutine elj_long(evdw) @@ -11255,7 +11993,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTACTS' real(kind=8),parameter :: accur=1.0d-10 - real(kind=8),dimension(3) :: gg + real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj !el local variables integer :: i,iint,j,k,itypi,itypi1,itypj real(kind=8) :: xi,yi,zi,xj,yj,zj,rij,sss,rrij,fac,eps0ij @@ -11287,8 +12025,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=e1+e2 evdw=evdw+(1.0d0-sss)*evdwij ! @@ -11345,7 +12083,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTACTS' real(kind=8),parameter :: accur=1.0d-10 - real(kind=8),dimension(3) :: gg + real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj !el local variables integer :: i,iint,j,k,itypi,itypi1,itypj,num_conti real(kind=8) :: xi,yi,zi,xj,yj,zj,rij,sss,rrij,fac,eps0ij @@ -11380,8 +12118,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=e1+e2 evdw=evdw+sss*evdwij ! @@ -11434,7 +12172,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' - real(kind=8),dimension(3) :: gg + real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj logical :: scheck !el local variables integer :: i,iint,j,k,itypi,itypi1,itypj @@ -11468,8 +12206,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (sss.lt.1.0d0) then 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=e_augm+e1+e2 !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -11521,7 +12259,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' - real(kind=8),dimension(3) :: gg + real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj logical :: scheck !el local variables integer :: i,iint,j,k,itypi,itypi1,itypj @@ -11555,8 +12293,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (sss.gt.0.0d0) then 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=e_augm+e1+e2 !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -11676,16 +12414,16 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Calculate whole angle-dependent part of epsilon and contributions ! to its derivatives fac=(rrij*sigsq)**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*(1.0d0-sss) 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) !d write (iout,'(2(a3,i3,2x),15(0pf7.3))') !d & restyp(itypi),i,restyp(itypj),j, !d & epsi,sigm,chi1,chi2,chip1,chip2, @@ -11796,16 +12534,16 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Calculate whole angle-dependent part of epsilon and contributions ! to its derivatives fac=(rrij*sigsq)**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*sss 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) !d write (iout,'(2(a3,i3,2x),15(0pf7.3))') !d & restyp(itypi),i,restyp(itypj),j, !d & epsi,sigm,chi1,chi2,chip1,chip2, @@ -11857,7 +12595,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift real(kind=8) :: sss,e1,e2,evdw,sss_grad real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& - dist_temp, dist_init + dist_temp, dist_init,aa,bb,fracinbuf,sslipi,sslipj,& + ssgradlipi,ssgradlipj + evdw=0.0D0 !cccc energy_dec=.false. @@ -11879,6 +12619,29 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) dzi=dc_norm(3,nres+i) @@ -11891,6 +12654,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & + 'evdw',i,j,evdwij,' ss' +! if (energy_dec) write (iout,*) & +! 'evdw',i,j,evdwij,' ss' + ELSE !el ind=ind+1 itypj=itype(j) if (itypj.eq.ntyp1) cycle @@ -11919,6 +12690,33 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj @@ -11983,8 +12781,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !--------------------------------------------------------------- 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 @@ -11993,8 +12791,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*(1.0d0-sss)*sss_ele_cut 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j,& epsi,sigm,chi1,chi2,chip1,chip2,& @@ -12023,6 +12821,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gg(3)=zj*fac ! Calculate angular part of the gradient. call sc_grad_scale(1.0d0-sss) + ENDIF !mask_dyn_ss endif enddo ! j enddo ! iint @@ -12056,7 +12855,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& - dist_temp, dist_init + dist_temp, dist_init,aa,bb,fracinbuf,sslipi,sslipj,& + ssgradlipi,ssgradlipj evdw=0.0D0 !cccc energy_dec=.false. ! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -12077,6 +12877,29 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) dzi=dc_norm(3,nres+i) @@ -12095,6 +12918,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & + 'evdw',i,j,evdwij,' ss' +! if (energy_dec) write (iout,*) & +! 'evdw',i,j,evdwij,' ss' + ELSE !el ind=ind+1 itypj=itype(j) if (itypj.eq.ntyp1) cycle @@ -12126,11 +12957,39 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj zj_safe=zj subchap=0 + do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 @@ -12191,8 +13050,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !--------------------------------------------------------------- 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 @@ -12201,8 +13060,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*sss*sss_ele_cut 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j,& epsi,sigm,chi1,chi2,chip1,chip2,& @@ -12233,6 +13092,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Calculate angular part of the gradient. call sc_grad_scale(sss) endif + ENDIF !mask_dyn_ss enddo ! j enddo ! iint enddo ! i @@ -12333,8 +13193,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !--------------------------------------------------------------- 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -12343,8 +13203,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' evdwij=evdwij*eps2rt*eps3rt evdw=evdw+(evdwij+e_augm)*(1.0d0-sss) 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j,& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),& @@ -12462,8 +13322,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !--------------------------------------------------------------- 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -12472,8 +13332,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' evdwij=evdwij*eps2rt*eps3rt evdw=evdw+(evdwij+e_augm)*sss 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j,& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),& @@ -12584,7 +13444,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #ifdef TIMING time01=MPI_Wtime() #endif +! print *, "before set matrices" call set_matrices +! print *,"after set martices" #ifdef TIMING time_mat=time_mat+MPI_Wtime()-time01 #endif @@ -14880,7 +15742,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! include 'COMMON.SCCOR' ! !el local variables - integer :: i,j,intertyp + integer :: i,j,intertyp,k ! Initialize Cartesian-coordinate gradient ! ! if (.not.allocated(gradx)) allocate(gradx(3,nres,2)) !(3,maxres,2) @@ -14922,7 +15784,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! allocate(gloc_sc(3,nres,10)) !(3,0:maxres2,10)maxres2=2*maxres !elwrite(iout,*) "icg",icg - do i=1,nres + do i=-1,nres do j=1,3 gvdwx(j,i)=0.0D0 gradx_scp(j,i)=0.0D0 @@ -14954,11 +15816,40 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gradx(j,i,icg)=0.0d0 gscloc(j,i)=0.0d0 gsclocx(j,i)=0.0d0 + gliptran(j,i)=0.0d0 + gshieldx(j,i)=0.0d0 + gshieldc(j,i)=0.0d0 + gshieldc_loc(j,i)=0.0d0 + gshieldx_ec(j,i)=0.0d0 + gshieldc_ec(j,i)=0.0d0 + gshieldc_loc_ec(j,i)=0.0d0 + gshieldx_t3(j,i)=0.0d0 + gshieldc_t3(j,i)=0.0d0 + gshieldc_loc_t3(j,i)=0.0d0 + gshieldx_t4(j,i)=0.0d0 + gshieldc_t4(j,i)=0.0d0 + gshieldc_loc_t4(j,i)=0.0d0 + gshieldx_ll(j,i)=0.0d0 + gshieldc_ll(j,i)=0.0d0 + gshieldc_loc_ll(j,i)=0.0d0 + do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo enddo enddo + do i=1,nres + do j=1,maxcontsshi + shield_list(j,i)=0 + do k=1,3 +!C print *,i,j,k + grad_shield_side(k,j,i)=0.0d0 + grad_shield_loc(k,j,i)=0.0d0 + enddo + enddo + ishield_list(i)=0 + enddo + ! ! Initialize the gradient of local energy terms. ! @@ -16418,9 +17309,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ljXs=sig-sig0ij ljA=eps1*eps2rt**2*eps3rt**2 - ljB=ljA*bb(itypi,itypj) - ljA=ljA*aa(itypi,itypj) - ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) + ljB=ljA*bb_aq(itypi,itypj) + ljA=ljA*aa_aq(itypi,itypj) + ljxm=ljXs+(-2.0D0*aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) ssXs=d0cm deltat1=1.0d0-om1 @@ -16454,7 +17345,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Stop and plot energy and derivative as a function of distance if (checkstop) then ssm=ssC-0.25D0*ssB*ssB/ssA - ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) + ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj) if (ssm.lt.ljm .and. & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then nicheck=1000 @@ -16479,8 +17370,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' havebond=.false. ljd=rij-ljXs fac=(1.0D0/ljd)**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) eij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=eij*eps3rt eps3der=eij*eps2rt @@ -16545,8 +17436,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3) else havebond=.false. - ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) - d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB + ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj) + d_ljm(1)=-0.5D0*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)*ljB d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt) d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- & alf12/eps3rt) @@ -16878,6 +17769,292 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' return end subroutine dyn_set_nss +! Lipid transfer energy function + subroutine Eliptransfer(eliptran) +!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 + real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip + integer :: i + eliptran=0.0 + print *, "I am in eliptran" + do i=ilip_start,ilip_end +!C do i=1,1 + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1).or.(i.eq.nres))& + cycle + + positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize)) + if (positi.le.0.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 + +!C print *,"doing sccale for lower part" +!C print *,i,sslip,fracinbuf,ssgradlip + 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 + if (energy_dec) write(iout,*) i,"eliptran=",eliptran,positi,sslip + enddo +! here starts the side chain transfer + do i=ilip_start,ilip_end + 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 + if (energy_dec) write(iout,*) i,"eliptran=",eliptran + enddo + return + end subroutine Eliptransfer +!-------------------------------------------------------------------------------- +!C first for shielding is setting of function of side-chains + + subroutine set_shield_fac2 + real(kind=8) :: div77_81=0.974996043d0, & + div4_81=0.2222222222d0 + real (kind=8) :: dist_pep_side,dist_side_calf,dist_pept_group, & + scale_fac_dist,fac_help_scale,VofOverlap,VolumeTotal,costhet,& + short,long,sinthet,costhet_fac,sh_frac_dist,rkprim,cosphi, & + sinphi,cosphi_fac,pep_side0pept_group,cosalfa,fac_alfa_sin +!C the vector between center of side_chain and peptide group + real(kind=8),dimension(3) :: pep_side_long,side_calf, & + pept_group,costhet_grad,cosphi_grad_long, & + cosphi_grad_loc,pep_side_norm,side_calf_norm, & + sh_frac_dist_grad,pep_side + integer i,j,k +!C write(2,*) "ivec",ivec_start,ivec_end + do i=1,nres + fac_shield(i)=0.0d0 + do j=1,3 + grad_shield(j,i)=0.0d0 + enddo + enddo + do i=ivec_start,ivec_end +!C do i=1,nres-1 +!C if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle + ishield_list(i)=0 + if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle +!Cif there two consequtive dummy atoms there is no peptide group between them +!C the line below has to be changed for FGPROC>1 + VolumeTotal=0.0 + do k=1,nres + if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle + dist_pep_side=0.0 + dist_side_calf=0.0 + do j=1,3 +!C first lets set vector conecting the ithe side-chain with kth side-chain + pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 +!C pep_side(j)=2.0d0 +!C and vector conecting the side-chain with its proper calfa + side_calf(j)=c(j,k+nres)-c(j,k) +!C side_calf(j)=2.0d0 + pept_group(j)=c(j,i)-c(j,i+1) +!C lets have their lenght + dist_pep_side=pep_side(j)**2+dist_pep_side + dist_side_calf=dist_side_calf+side_calf(j)**2 + dist_pept_group=dist_pept_group+pept_group(j)**2 + enddo + dist_pep_side=sqrt(dist_pep_side) + dist_pept_group=sqrt(dist_pept_group) + dist_side_calf=sqrt(dist_side_calf) + do j=1,3 + pep_side_norm(j)=pep_side(j)/dist_pep_side + side_calf_norm(j)=dist_side_calf + enddo +!C now sscale fraction + sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield +!C print *,buff_shield,"buff" +!C now sscale + if (sh_frac_dist.le.0.0) cycle +!C print *,ishield_list(i),i +!C If we reach here it means that this side chain reaches the shielding sphere +!C Lets add him to the list for gradient + ishield_list(i)=ishield_list(i)+1 +!C ishield_list is a list of non 0 side-chain that contribute to factor gradient +!C this list is essential otherwise problem would be O3 + shield_list(ishield_list(i),i)=k +!C Lets have the sscale value + if (sh_frac_dist.gt.1.0) then + scale_fac_dist=1.0d0 + do j=1,3 + sh_frac_dist_grad(j)=0.0d0 + enddo + else + scale_fac_dist=-sh_frac_dist*sh_frac_dist & + *(2.0d0*sh_frac_dist-3.0d0) + fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2) & + /dist_pep_side/buff_shield*0.5d0 + do j=1,3 + sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) +!C sh_frac_dist_grad(j)=0.0d0 +!C scale_fac_dist=1.0d0 +!C print *,"jestem",scale_fac_dist,fac_help_scale, +!C & sh_frac_dist_grad(j) + enddo + endif +!C this is what is now we have the distance scaling now volume... + short=short_r_sidechain(itype(k)) + long=long_r_sidechain(itype(k)) + costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2) + sinthet=short/dist_pep_side*costhet +!C now costhet_grad +!C costhet=0.6d0 +!C sinthet=0.8 + costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4 +!C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet +!C & -short/dist_pep_side**2/costhet) +!C costhet_fac=0.0d0 + do j=1,3 + costhet_grad(j)=costhet_fac*pep_side(j) + enddo +!C remember for the final gradient multiply costhet_grad(j) +!C for side_chain by factor -2 ! +!C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 +!C pep_side0pept_group is vector multiplication + pep_side0pept_group=0.0d0 + do j=1,3 + pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) + enddo + cosalfa=(pep_side0pept_group/ & + (dist_pep_side*dist_side_calf)) + fac_alfa_sin=1.0d0-cosalfa**2 + fac_alfa_sin=dsqrt(fac_alfa_sin) + rkprim=fac_alfa_sin*(long-short)+short +!C rkprim=short + +!C now costhet_grad + cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2) +!C cosphi=0.6 + cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4 + sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/ & + dist_pep_side**2) +!C sinphi=0.8 + do j=1,3 + cosphi_grad_long(j)=cosphi_fac*pep_side(j) & + +cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) & + *(long-short)/fac_alfa_sin*cosalfa/ & + ((dist_pep_side*dist_side_calf))* & + ((side_calf(j))-cosalfa* & + ((pep_side(j)/dist_pep_side)*dist_side_calf)) +!C cosphi_grad_long(j)=0.0d0 + cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) & + *(long-short)/fac_alfa_sin*cosalfa & + /((dist_pep_side*dist_side_calf))* & + (pep_side(j)- & + cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) +!C cosphi_grad_loc(j)=0.0d0 + enddo +!C print *,sinphi,sinthet + VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) & + & /VSolvSphere_div +!C & *wshield +!C now the gradient... + do j=1,3 + grad_shield(j,i)=grad_shield(j,i) & +!C gradient po skalowaniu + +(sh_frac_dist_grad(j)*VofOverlap & +!C gradient po costhet + +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0* & + (1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( & + sinphi/sinthet*costhet*costhet_grad(j) & + +sinthet/sinphi*cosphi*cosphi_grad_long(j))) & + )*wshield +!C grad_shield_side is Cbeta sidechain gradient + grad_shield_side(j,ishield_list(i),i)=& + (sh_frac_dist_grad(j)*-2.0d0& + *VofOverlap& + -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*& + (1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(& + sinphi/sinthet*costhet*costhet_grad(j)& + +sinthet/sinphi*cosphi*cosphi_grad_long(j))) & + )*wshield + + grad_shield_loc(j,ishield_list(i),i)= & + scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*& + (1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(& + sinthet/sinphi*cosphi*cosphi_grad_loc(j)& + ))& + *wshield + enddo + VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist + enddo + fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield) + + write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i) + enddo + return + end subroutine set_shield_fac2 + !----------------------------------------------------------------------------- #ifdef WHAM subroutine read_ssHist @@ -16984,9 +18161,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(grij_hb_cont(3,maxconts,nres)) !(3,maxconts,maxres) allocate(facont_hb(maxconts,nres)) + allocate(ees0p(maxconts,nres)) allocate(ees0m(maxconts,nres)) allocate(d_cont(maxconts,nres)) + allocate(ees0plist(maxconts,nres)) + !(maxconts,maxres) allocate(num_cont_hb(nres)) !(maxres) @@ -17079,38 +18259,60 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !(6,maxdim) allocate(dxds(6,nres)) !(6,maxres) - allocate(gradx(3,nres,0:2)) - allocate(gradc(3,nres,0:2)) + allocate(gradx(3,-1:nres,0:2)) + allocate(gradc(3,-1:nres,0:2)) !(3,maxres,2) - allocate(gvdwx(3,nres)) - allocate(gvdwc(3,nres)) - allocate(gelc(3,nres)) - allocate(gelc_long(3,nres)) - allocate(gvdwpp(3,nres)) - allocate(gvdwc_scpp(3,nres)) - allocate(gradx_scp(3,nres)) - allocate(gvdwc_scp(3,nres)) - allocate(ghpbx(3,nres)) - allocate(ghpbc(3,nres)) - allocate(gradcorr(3,nres)) - allocate(gradcorr_long(3,nres)) - allocate(gradcorr5_long(3,nres)) - allocate(gradcorr6_long(3,nres)) - allocate(gcorr6_turn_long(3,nres)) - allocate(gradxorr(3,nres)) - allocate(gradcorr5(3,nres)) - allocate(gradcorr6(3,nres)) + allocate(gvdwx(3,-1:nres)) + allocate(gvdwc(3,-1:nres)) + allocate(gelc(3,-1:nres)) + allocate(gelc_long(3,-1:nres)) + allocate(gvdwpp(3,-1:nres)) + allocate(gvdwc_scpp(3,-1:nres)) + allocate(gradx_scp(3,-1:nres)) + allocate(gvdwc_scp(3,-1:nres)) + allocate(ghpbx(3,-1:nres)) + allocate(ghpbc(3,-1:nres)) + allocate(gradcorr(3,-1:nres)) + allocate(gradcorr_long(3,-1:nres)) + allocate(gradcorr5_long(3,-1:nres)) + allocate(gradcorr6_long(3,-1:nres)) + allocate(gcorr6_turn_long(3,-1:nres)) + allocate(gradxorr(3,-1:nres)) + allocate(gradcorr5(3,-1:nres)) + allocate(gradcorr6(3,-1:nres)) + allocate(gliptran(3,-1:nres)) + allocate(gliptranc(3,-1:nres)) + allocate(gliptranx(3,-1:nres)) + allocate(gshieldx(3,-1:nres)) + allocate(gshieldc(3,-1:nres)) + allocate(gshieldc_loc(3,-1:nres)) + allocate(gshieldx_ec(3,-1:nres)) + allocate(gshieldc_ec(3,-1:nres)) + allocate(gshieldc_loc_ec(3,-1:nres)) + allocate(gshieldx_t3(3,-1:nres)) + allocate(gshieldc_t3(3,-1:nres)) + allocate(gshieldc_loc_t3(3,-1:nres)) + allocate(gshieldx_t4(3,-1:nres)) + allocate(gshieldc_t4(3,-1:nres)) + allocate(gshieldc_loc_t4(3,-1:nres)) + allocate(gshieldx_ll(3,-1:nres)) + allocate(gshieldc_ll(3,-1:nres)) + allocate(gshieldc_loc_ll(3,-1:nres)) + allocate(grad_shield(3,-1:nres)) !(3,maxres) + allocate(grad_shield_side(3,50,nres)) + allocate(grad_shield_loc(3,50,nres)) +! grad for shielding surroing allocate(gloc(0:maxvar,0:2)) allocate(gloc_x(0:maxvar,2)) !(maxvar,2) - allocate(gel_loc(3,nres)) - allocate(gel_loc_long(3,nres)) - allocate(gcorr3_turn(3,nres)) - allocate(gcorr4_turn(3,nres)) - allocate(gcorr6_turn(3,nres)) - allocate(gradb(3,nres)) - allocate(gradbx(3,nres)) + allocate(gel_loc(3,-1:nres)) + allocate(gel_loc_long(3,-1:nres)) + allocate(gcorr3_turn(3,-1:nres)) + allocate(gcorr4_turn(3,-1:nres)) + allocate(gcorr6_turn(3,-1:nres)) + allocate(gradb(3,-1:nres)) + allocate(gradbx(3,-1:nres)) !(3,maxres) allocate(gel_loc_loc(maxvar)) allocate(gel_loc_turn3(maxvar)) @@ -17120,19 +18322,19 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(g_corr5_loc(maxvar)) allocate(g_corr6_loc(maxvar)) !(maxvar) - allocate(gsccorc(3,nres)) - allocate(gsccorx(3,nres)) + allocate(gsccorc(3,-1:nres)) + allocate(gsccorx(3,-1:nres)) !(3,maxres) - allocate(gsccor_loc(nres)) + allocate(gsccor_loc(-1:nres)) !(maxres) - allocate(dtheta(3,2,nres)) + allocate(dtheta(3,2,-1:nres)) !(3,2,maxres) - allocate(gscloc(3,nres)) - allocate(gsclocx(3,nres)) + allocate(gscloc(3,-1:nres)) + allocate(gsclocx(3,-1:nres)) !(3,maxres) - allocate(dphi(3,3,nres)) - allocate(dalpha(3,3,nres)) - allocate(domega(3,3,nres)) + allocate(dphi(3,3,-1:nres)) + allocate(dalpha(3,3,-1:nres)) + allocate(domega(3,3,-1:nres)) !(3,3,maxres) ! common /deriv_scloc/ allocate(dXX_C1tab(3,nres)) @@ -17170,11 +18372,11 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !---------------------- ! common.MD ! common /mdgrad/ - allocate(gcart(3,0:nres)) - allocate(gxcart(3,0:nres)) + allocate(gcart(3,-1:nres)) + allocate(gxcart(3,-1:nres)) !(3,0:MAXRES) - allocate(gradcag(3,nres)) - allocate(gradxag(3,nres)) + allocate(gradcag(3,-1:nres)) + allocate(gradxag(3,-1:nres)) !(3,MAXRES) ! common /back_constr/ !el in energy:Econstr_back allocate((:),allocatable :: utheta,ugamma,uscdiff !(maxfrag_back) @@ -17220,7 +18422,10 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(idssb(nss),jdssb(nss)) !(maxdim) endif + allocate(ishield_list(nres)) + allocate(shield_list(50,nres)) allocate(dyn_ss_mask(nres)) + allocate(fac_shield(nres)) !(maxres) dyn_ss_mask(:)=.false. !---------------------- diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 243c8b6..7964cdd 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -730,6 +730,9 @@ call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) + call reada(weightcard,'WSHIELD',wshield,0.05D0) + call reada(weightcard,'LIPSCALE',lipscale,1.0D0) + call reada(weightcard,'WLT',wliptran,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 451731d..0dce2e1 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -743,7 +743,7 @@ real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,& dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,& - res1 + res1,epsijlip integer :: ichir1,ichir2 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) @@ -782,7 +782,8 @@ allocate(isc(ntyp+1)) !(ntyp+1) allocate(restok(ntyp+1)) !(ntyp+1) allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) - + allocate(long_r_sidechain(ntyp)) + allocate(short_r_sidechain(ntyp)) dsc(:)=0.0d0 dsc_inv(:)=0.0d0 @@ -799,7 +800,7 @@ endif enddo #else - read (ibond,*) junk,vbldpDUM,vbldp0,akp,rjunk,mp,ip,pstok + read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok do i=1,ntyp read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),& j=1,nbondterm(i)),msc(i),isc(i),restok(i) @@ -841,6 +842,17 @@ theta0(:)=0.0D0 sig0(:)=0.0D0 sigc0(:)=0.0D0 + allocate(liptranene(ntyp)) +!C reading lipid parameters + write (iout,*) "iliptranpar",iliptranpar + call flush(iout) + read(iliptranpar,*) pepliptran + print *,pepliptran + do i=1,ntyp + read(iliptranpar,*) liptranene(i) + print *,liptranene(i) + enddo + close(iliptranpar) #ifdef CRYST_THETA ! @@ -1418,8 +1430,8 @@ do i=-ntyp,-1 itortyp(i)=-itortyp(-i) enddo -! itortyp(ntyp1)=ntortyp -! itortyp(-ntyp1)=-ntortyp + itortyp(ntyp1)=ntortyp + itortyp(-ntyp1)=-ntortyp do iblock=1,2 write (iout,*) 'ntortyp',ntortyp do i=0,ntortyp-1 @@ -1916,7 +1928,7 @@ allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2) allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp) allocate(chip(ntyp1),alp(ntyp1)) !(ntyp) - + allocate(epslip(ntyp,ntyp)) augm(:,:)=0.0D0 chip(:)=0.0D0 alp(:)=0.0D0 @@ -1980,6 +1992,10 @@ read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp) read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp) read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp) + do i=1,ntyp + read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp) + enddo + ! For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -2023,10 +2039,13 @@ ! Calculate the "working" parameters of SC interactions. !el from module energy - COMMON.INTERACT------- - allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp) allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) - aa(:,:)=0.0D0 - bb(:,:)=0.0D0 + aa_aq(:,:)=0.0D0 + bb_aq(:,:)=0.0D0 + aa_lip(:,:)=0.0D0 + bb_lip(:,:)=0.0D0 chi(:,:)=0.0D0 sigma(:,:)=0.0D0 r0(:,:)=0.0D0 @@ -2036,6 +2055,7 @@ 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 @@ -2064,10 +2084,18 @@ 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) + epsijlip=epslip(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) +!C write(iout,*) aa_lip if (ipot.gt.2) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 @@ -2100,7 +2128,7 @@ endif if (lprint) then write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') & - restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),& + restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),& sigma(i,j),r0(i,j),chi(i,j),chi(j,i) endif enddo @@ -2523,11 +2551,11 @@ e2(3)=0.0d0 endif do j=1,3 - c(j,nres)=c(j,nres-1)-3.8d0*e2(j) + c(j,nres)=c(j,nres-1)-1.9d0*e2(j) enddo else do j=1,3 - dcj=c(j,nres-2)-c(j,nres-3) + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo @@ -2592,11 +2620,11 @@ e2(3)=0.0d0 endif do j=1,3 - c(j,1)=c(j,2)-3.8d0*e2(j) + c(j,1)=c(j,2)-1.9d0*e2(j) enddo else do j=1,3 - dcj=c(j,4)-c(j,3) + dcj=(c(j,4)-c(j,3))/2.0 c(j,1)=c(j,2)-dcj c(j,nres+1)=c(j,1) enddo @@ -2747,6 +2775,9 @@ ! enddiagnostic ! makes copy of chains write (iout,*) "symetr", symetr + do j=1,3 + dc(j,0)=c(j,1) + enddo if (symetr.gt.1) then call permut(symetr) @@ -2861,7 +2892,7 @@ character(len=640) :: controlcard real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,& - + integer i nglob_csa=0 eglob_csa=1d99 @@ -2900,6 +2931,14 @@ timem=timlim modecalc=0 call reada(controlcard,"T_BATH",t_bath,300.0d0) +!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 +!C solvent accesible area +!C 1 the newly introduced function +!C 2 reseved for further possible developement + call readi(controlcard,'SHIELD',shield_mode,0) +!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "shield_mode",shield_mode !C Varibles set size of box call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) @@ -2907,6 +2946,25 @@ ! CUTOFFF ON ELECTROSTATICS call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0) + write(iout,*) "R_CUT_ELE=",r_cut_ele +! Lipidic parameters + 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 + 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 !lipthick.gt.0 + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot + write (iout,*) "SHIELD MODE",shield_mode !C------------------------- minim=(index(controlcard,'MINIMIZE').gt.0) @@ -3008,6 +3066,25 @@ if(me.eq.king.or..not.out1file) & write (iout,'(2a)') diagmeth(kdiag),& ' routine used to diagonalize matrices.' + if (shield_mode.gt.0) then + pi=3.141592d0 +!C VSolvSphere the volume of solving sphere +!C print *,pi,"pi" +!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +!C there will be no distinction between proline peptide group and normal peptide +!C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +!C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) + short_r_sidechain(i)=sigma0(i) + write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& + sigma0(i) + enddo + buff_shield=1.0d0 + endif return end subroutine read_control !----------------------------------------------------------------------------- @@ -3666,6 +3743,8 @@ ! print *,"Processor",myrank," opened file IELEP" call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') ! print *,"Processor",myrank," opened file ISIDEP" ! print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -3691,6 +3770,8 @@ open (ielep,file=elename,status='old') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old') + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -3721,6 +3802,8 @@ open (ielep,file=elename,status='old',readonly) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') diff --git a/source/unres/minim.f90 b/source/unres/minim.f90 index d556363..269b904 100644 --- a/source/unres/minim.f90 +++ b/source/unres/minim.f90 @@ -4986,16 +4986,17 @@ !--------------------------------------------------------------- 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_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij + if (wliptran.gt.0.0) print *,"WARNING eps_aq used!" 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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') & !d restyp(itypi),i,restyp(itypj),j, & !d epsi,sigm,chi1,chi2,chip1,chip2, & -- 1.7.9.5