.0572167103 -.0468608825 .0151048455 .0084963678
.0278930397 .0076922911 .1033536738 -.0098256036
.0611385674 .0448303346 .0861379100 .0861379100
-C From here are lipid
- 5.6053537261 6.2200032154 6.2159898496 6.4386968963
- 6.2098101231 5.9596374798 5.4310222662 4.8790085257
- 5.2638423822 5.4961954124 4.2928124162 4.3582250540
- 4.2667755921 3.6813845958 3.5405429647 3.7026985912
- 4.7177596156 3.3811992227 3.7936085458 4.5157490585
- 6.2200032154 6.2159898496 5.2638423822 5.2638423822
-
- 6.6296746680 6.6717260508 6.8168876945 6.8378960152
- 6.3512039986 6.1425192423 5.5200206406 5.4858936024
- 4.9382573262 4.2795129166 4.0534857018 4.2086428525
- 3.5388086279 3.4209730327 2.7995744193 4.8164268553
- 3.8790755071 3.5878262177 4.6457432555 6.6296746680
- 6.6717260508 5.4858936024 5.4858936024
-
- 6.6424306340 6.9715947250 6.9241225787 6.5291325353
- 6.1394821460 5.4415971840 5.2914044780 4.7881880526
- 4.1302408718 3.9405275117 3.9022731880 3.6078928145
- 3.0809713206 3.0182595342 4.7935507272 3.8711928134
- 3.6826168780 4.4872661030 6.6717260508 6.6424306340
- 5.2914044780 5.2914044780
-
- 7.3271707934 7.2970176793 7.0494948235 6.2565687662
- 5.8133715113 5.8845382769 5.3463553241 4.7146928716
- 4.4864793359 4.3092352779 3.6928361281 4.0100566871
- 3.6584608489 4.5901676001 4.4681844010 4.6334334282
- 4.9023945082 6.8168876945 6.9715947250 5.8845382769
- 5.8845382769
-
- 7.0106342493 6.8341607204 6.2047876082 5.5802122247
- 5.6716968726 5.2011344932 4.1966522349 4.2383931955
- 3.8954073453 3.5679086118 3.1136348127 2.7622561241
- 4.5092053374 4.0238412606 3.9799853578 4.8087180377
- 6.8378960152 6.9241225787 5.6716968726 5.6716968726
-
- 6.5353606708 5.8844482033 5.1100976261 5.5515360893
- 4.9587559238 4.2030790774 4.0351600709 3.6506289806
- 3.4640416936 3.1611022087 2.7976192835 4.0422832763
- 3.4415952740 3.8620694983 4.5934246009 6.3512039986
- 6.5291325353 5.5515360893 5.5515360893
-
- 5.2828811325 4.8298466315 4.7575327777 4.7249779042
- 3.5656384285 3.5613366567 3.9273045564 3.7664779118
- 3.5861699761 3.5523191627 4.6508672463 4.2424308260
- 4.0987727783 4.5079669929 6.1425192423 6.1394821460
- 4.7575327777 4.7575327777
-
- 4.2222645754 4.1627149417 3.9786372270 3.1623319201
- 2.9379696734 3.1365985188 2.8025851476 2.6576343487
- 2.8733395837 3.8785891655 3.4443068252 3.4405132673
- 4.0915305763 5.5200206406 5.4415971840 4.1627149417
- 4.1627149417
-
- 4.3860654351 3.6520014167 3.0264937881 2.5680579493
- 2.5001857370 2.1393094782 1.4707077389 1.7373181916
- 3.1045342752 2.0152329038 1.8196209136 3.6395852546
- 5.4858936024 5.2914044780 4.3860654351 4.3860654351
-
- 2.5016557935 2.7879606493 2.1424756816 1.0741544747
- 1.3236438850 -.2712356678 .7613312181 2.8627348285
- 1.5002127649 -.2392577636 3.5381128985 4.9382573262
- 4.7881880526 3.6520014167 3.6520014167
-
- 2.2903847447 2.2466771550 1.5053607856 1.4031797188
- .9676209687 1.3075541842 2.8595118604 1.9833053538
- -.0085773647 2.9557553666 4.2795129166 4.1302408718
- 3.0264937881 3.0264937881
-
- 1.2800038243 .7689806100 .9119606366 -.0744371793
- .3669965427 2.4288563200 1.4975482174 -.2134101897
- 2.9415021479 4.0534857018 3.9405275117 2.5680579493
- 2.5680579493
-
- -.6792428859 .4532383239 -.7590387660 -.3617034846
- 1.6803275058 .6775209988 -.5354837468 2.6208591363
- 4.2086428525 3.9022731880 2.5001857370 2.5001857370
-
- .4136122500 -.1637152873 .2212564728 2.1287126020
- .3896388275 -.2354089650 2.3255326077 3.5388086279
- 3.6078928145 2.1393094782 2.1393094782
-
- -3.5193531074 -2.0157307929 1.8355926160 2.7489174124
- 2.0098391558 1.7975718667 3.4209730327 3.0809713206
- 1.4707077389 1.4707077389
-
- -1.3532688232 2.5504275249 2.8202873790 1.6253609671
- 1.8625091247 2.7995744193 3.0182595342 1.7373181916
- 1.7373181916
-
- 3.7292778697 2.2944436481 -.0303565255 3.1115776177
- 4.8164268553 4.7935507272 3.1045342752 3.1045342752
-
- -.0827362961 -1.6043113182 2.4439837435 3.8790755071
- 3.8711928134 2.0152329038 2.0152329038
-
- -4.3897783280 2.3664634533 3.5878262177 3.6826168780
- 1.8196209136 1.8196209136
-
- 4.1927969260 4.6457432555 4.4872661030 3.6395852546
- 3.6395852546
-
- 6.6296746680 6.6717260508 5.4858936024 5.4858936024
-
- 6.6424306340 5.2914044780 5.2914044780
-
- 4.3860654351 4.3860654351
-
- 4.3860654351
-
+ 2.252000 2.758289 2.828747 2.573807
+ 2.573807 2.314337 3.329045 3.123838
+ 1.679138 1.245463 2.417517 2.138542
+ 2.902927 2.675567 2.916240 2.690005
+ 2.934561 3.360153 2.904943 2.252400
+ 2.822849 3.512073 2.021838 2.021838
+
+ 3.378400 3.464698 3.152444 3.152444
+ 2.834641 4.077473 3.826132 2.056637
+ 1.525464 2.961017 2.619323 3.555555
+ 3.277081 3.571861 3.294765 3.594301
+ 4.115573 3.558025 2.758779 3.457474
+ 4.301648 2.476383 2.476383
+
+ 3.553200 3.232970 3.232970 2.907049
+ 4.181627 3.923866 2.109172 1.564431
+ 3.036654 2.686231 3.646378 3.360790
+ 3.663100 3.378926 3.686114 4.220702
+ 3.648911 2.829249 3.545792 4.411529
+ 2.539639 2.539639
+
+ 2.941600 2.941600 2.645052 3.804760
+ 3.570230 1.919084 1.423437 2.762977
+ 2.444136 3.317750 3.057901 3.332965
+ 3.074402 3.353905 3.840313 3.320055
+ 2.574264 3.226230 4.013942 2.310756
+ 2.310756
+
+ 2.941600 2.645052 3.804760 3.570230
+ 1.919084 1.423437 2.762977 2.444136
+ 3.317750 3.057901 3.332965 3.074402
+ 3.353905 3.840313 3.320055 2.574264
+ 3.226230 4.013942 2.310756 2.310756
+
+ 2.378400 3.421196 3.210309 1.725618
+ 1.279938 2.484436 2.197739 2.983282
+ 2.749629 2.996964 2.764467 3.015792
+ 3.453165 2.985354 2.314748 2.900988
+ 3.609290 2.077805 2.077805
+
+ 4.921200 4.617850 2.482205 1.841120
+ 3.573723 3.161325 4.291286 3.955188
+ 4.310965 3.976531 4.338049 4.967185
+ 4.294266 3.329637 4.172909 5.191762
+ 2.988806 2.988806
+
+ 4.333200 2.329199 1.727631 3.353434
+ 2.966456 4.026765 3.711385 4.045231
+ 3.731412 4.070646 4.661001 4.029562
+ 3.124393 3.915686 4.871735 2.804572
+ 2.804572
+
+ 1.252000 0.9286429 1.802551 1.594541
+ 2.164482 1.994958 2.174409 2.005723
+ 2.188069 2.505400 2.165986 1.679436
+ 2.104775 2.618674 1.507524 1.507524
+
+ 0.6888000 1.337002 1.182715 1.605456
+ 1.479715 1.612819 1.487700 1.622951
+ 1.858324 1.606571 1.245684 1.561169
+ 1.942343 1.118173 1.118173
+
+ 2.595200 2.295720 3.116286 2.872216
+ 3.130577 2.887715 3.150245 3.607117
+ 3.118450 2.417947 3.030323 3.770203
+ 2.170439 2.170439
+
+ 2.030800 2.756675 2.540769 2.769316
+ 2.554480 2.786715 3.190865 2.758589
+ 2.138922 2.680631 3.335131 1.919976
+ 1.919976
+
+ 3.742000 3.448923 3.759161 3.467535
+ 3.782778 4.331384 3.744599 2.903442
+ 3.638776 4.527216 2.606238 2.606238
+
+ 3.178800 3.464740 3.195954 3.486507
+ 3.992146 3.451319 2.676042 3.353784
+ 4.172640 2.402115 2.402115
+
+ -61.60000 -35.26000 3.800125 48.12500
+ 35.00000 2.916757 3.655463 4.547977
+ 2.618191 2.618191
+
+ -23.62500 3.505321 39.90000 28.44000
+ 2.690483 3.371882 4.195157 2.415078
+ 2.415078
+
+ 3.824000 4.378585 3.785405 2.935082
+ 3.678429 4.576550 2.634640 2.634640
+
+ -1.447000 -28.07000 3.360749 4.211903
+ 5.240276 3.016734 3.016734
+
+ -76.82000 2.905459 3.641304 4.530360
+ 2.608049 2.608049
+
+ 2.252800 2.823350 3.512697 2.022197
+ 2.022197
+
+ 3.538400 4.402332 2.534345 2.534345
+
+ 5.477200 3.153127 3.153127
+
+ 1.815200 1.815200
+
+ 1.815200
+
& "Time step reduced to",d_time,
& " because of too large initial acceleration."
endif
- if(me.eq.king.or..not.out1file)then
- write(iout,*) "Potential energy and its components"
- call enerprint(potEcomp)
+C if(me.eq.king.or..not.out1file)then
+C write(iout,*) "Potential energy and its components"
+C call enerprint(potEcomp)
c write(iout,*) (potEcomp(i),i=0,n_ene)
- endif
+C endif
potE=potEcomp(0)-potEcomp(20)
totE=EK+potE
itime=0
c write (iout,*) 'theta=', theta(i-1)
enddo
#else
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ iti = itortyp(itype(i-2))
+ else
+ iti=ntortyp+1
+ endif
+c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ iti1 = itortyp(itype(i-1))
+ else
+ iti1=ntortyp+1
+ endif
b1(1,i-2)=b(3,iti)
b1(2,i-2)=b(5,iti)
b2(1,i-2)=b(2,iti)
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
+C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
c write (iout,*) 'mu ',mu(:,i-2),i-2
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
& +a33*muij(4)
c write (iout,*) 'i',i,' j',j,itype(i),itype(j),
c & ' eel_loc_ij',eel_loc_ij
-c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
+C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
C Calculate patrial derivative for theta angle
#ifdef NEWCORR
geel_loc_ij=a22*gmuij1(1)
C now we start reading lipid
do i=1,ntyp
read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
- print *,"WARNING!!"
- do j=1,ntyp
- epslip(i,j)=epslip(i,j)+0.05d0
- enddo
+C print *,"WARNING!!"
+C do j=1,ntyp
+C epslip(i,j)=epslip(i,j)+0.05d0
+C enddo
enddo
+ write(iout,*) epslip(1,1),"OK?"
C For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do i=1,ntyp
bordliptop=(boxzsize+lipthick)/2.0
bordlipbot=bordliptop-lipthick
C endif
- if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0))
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
& write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
buflipbot=bordlipbot+lipbufthick
bufliptop=bordliptop-lipbufthick
& app_all(2,2,max_parm),bpp_all(2,2,max_parm),
& ael6_all(2,2,max_parm),ael3_all(2,2,max_parm),
& aad_all(ntyp,2,max_parm),bad_all(ntyp,2,max_parm),
- & aa_all(ntyp,ntyp,max_parm),bb_all(ntyp,ntyp,max_parm),
+ & aa_aq_all(ntyp,ntyp,max_parm),bb_aq_all(ntyp,ntyp,max_parm),
+ & aa_lip_all(ntyp,ntyp,max_parm),bb_lip_all(ntyp,ntyp,max_parm),
& augm_all(ntyp,ntyp,max_parm),eps_all(ntyp,ntyp,max_parm),
+ & epslip_all(ntyp,ntyp,max_parm),
& sigma_all(ntyp,ntyp,max_parm),r0_all(ntyp,ntyp,max_parm),
& chi_all(ntyp,ntyp,max_parm),chip_all(ntyp,max_parm),
& alp_all(ntyp,max_parm),ebr_all(max_parm),d0cm_all(max_parm),
& v0_all,v1_all,v2_all,vlor1_all,vlor2_all,vlor3_all,v1c_all,
& v1s_all,v2c_all,v2s_all,b1_all,b2_all,cc_all,dd_all,ee_all,
& ctilde_all,dtilde_all,b1tilde_all,app_all,bpp_all,ael6_all,
- & ael3_all,aad_all,bad_all,aa_all,bb_all,augm_all,
+ & ael3_all,aad_all,bad_all,aa_aq_all,bb_aq_all,augm_all,
+ & aa_lip_all,bb_lip_all,epslip_all,
& eps_all,sigma_all,r0_all,chi_all,chip_all,alp_all,ebr_all,
& d0cm_all,akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all,
& v1sccor_all,v2sccor_all,nbondterm_all,
&nend_sup,chain_length,tabperm(maxperm,maxsym),nperm,
& nstart_seq,ishift_pdb
double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,
- &sssgrad
- common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad
+ &sssgrad,
+ & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+ common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
+ & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
C General I/O units & files
integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam,
& itorp,itordp,ifourier,ielep,isidep,iscpp,isccor,icbase,
- & istat,ientin,ientout,isidep1,ibond,ihist,izsc,idistr
+ & istat,ientin,ientout,isidep1,ibond,ihist,izsc,idistr,
+ & iliptranpar
common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,
& irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,isccor,
& icbase,istat,ientin,ientout,isidep1,ibond,ihist,izsc,
- & idistr
+ & idistr,iliptranpar
character*256 outname,intname,pdbname,mol2name,statname,intinname,
& entname,restartname,prefix,scratchdir,sidepname,pdbfile,
& histname,zscname
& sidepname,pdbfile,histname,zscname
C Parameter files
character*256 bondname,thetname,rotname,torname,tordname,
- & fouriername,elename,sidename,scpname,sccorname,patname
+ & fouriername,elename,sidename,scpname,sccorname,patname,
+ & liptranname
common /parfiles/ thetname,rotname,torname,tordname,bondname,
- & fouriername,elename,sidename,scpname,sccorname,patname
+ & fouriername,elename,sidename,scpname,sccorname,patname,
+ & liptranname
character*3 pot
C-----------------------------------------------------------------------
C INP - main input file
endif
C write (iout,*) "Czy tu dochodze"
potE(iii+1,iparm)=energia(0)
- do k=1,21
+ do k=1,22
enetb(k,iii+1,iparm)=energia(k)
enddo
#ifdef DEBUG
C 21/5/07 Calculate local sicdechain correlation energy
C
call eback_sc_corr(esccor)
+
+ if (wliptran.gt.0) then
+ call Eliptransfer(eliptran)
+ endif
+
C
C 12/1/95 Multi-body terms
C
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+ & +wliptran*eliptran
#else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+ & +wliptran*eliptran
#endif
energia(0)=etot
energia(1)=evdw
energia(20)=edihcnstr
energia(21)=evdw_t
energia(24)=ethetacnstr
+ energia(22)=eliptran
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*fact(2)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
#else
do i=1,nct
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*fact(1)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
#endif
enddo
edihcnstr=energia(20)
estr=energia(18)
ethetacnstr=energia(24)
+ eliptran=energia(22)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
& wvdwpp,
& ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
& eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
& eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
- & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot
+ & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,
+ & eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond,
& ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
& eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
& eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
- & edihcnstr,ethetacnstr,ebr*nss,etot
+ & edihcnstr,ethetacnstr,ebr*nss,eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
ij=icant(itypi,itypj)
c ROZNICA z cluster
cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
rij=1.0D0/r_inv_ij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
ij=icant(itypi,itypj)
eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux
& /dabs(eps(itypi,itypj))
eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
endif
if (calc_grad) then
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),15(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
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)
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
C checking the distance
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
if (sss.le.0.0) cycle
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
+
call sc_angular
sigsq=1.0D0/sigsq
sig=sig0ij*dsqrt(sigsq)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0) then
+ if (bb.gt.0) then
evdw=evdw+evdwij*sss
else
evdw_t=evdw_t+evdwij*sss
c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
c & aux*e2/eps(itypi,itypj)
c if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
#ifdef DEBUG
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij+e_augm
else
evdw_t=evdw_t+evdwij+e_augm
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
enddo
+C write (iout,*) 'mumu',i,b1(1,iti),Ub2(1,i-2)
+
C Vectors and matrices dependent on a single virtual-bond dihedral.
call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2))
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
& +a33*muij(4)
c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
-c write (iout,'(a6,2i5,0pf7.3)')
-c & 'eelloc',i,j,eel_loc_ij
+C write (iout,'(a6,2i5,0pf7.3)')
+C & 'eelloc',i,j,eel_loc_ij
+C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
eel_loc=eel_loc+eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
enddo
endif
else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+5).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes suggested by Ana
+ & .or. itype(i+3).eq.ntyp1
+ & .or. itype(i+4).eq.ntyp1
+ & .or. itype(i+5).eq.ntyp1
+ & .or. itype(i).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1) goto 178
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Fourth-order contributions
gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
enddo
endif
+ 178 continue
endif
return
end
return
end
crc-------------------------------------------------
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+ subroutine Eliptransfer(eliptran)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+C this is done by Adasko
+C print *,"wchodze"
+C structure of box:
+C water
+C--bordliptop-- buffore starts
+C--bufliptop--- here true lipid starts
+C lipid
+C--buflipbot--- lipid ends buffore starts
+C--bordlipbot--buffore ends
+ eliptran=0.0
+ do i=1,nres
+C do i=1,1
+ if (itype(i).eq.ntyp1) cycle
+
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+C print *,i
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+ if ((positi.gt.bordlipbot)
+ &.and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+C print *, "doing sscalefor top part"
+C print *,i,sslip,fracinbuf,ssgradlip
+ else
+ eliptran=eliptran+pepliptran
+C print *,"I am in true lipid"
+ endif
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ endif
+ enddo
+C print *, "nic nie bylo w lipidzie?"
+C now multiply all by the peptide group transfer factor
+C eliptran=eliptran*pepliptran
+C now the same for side chains
+CV do i=1,1
+ do i=1,nres
+ if (itype(i).eq.ntyp1) cycle
+ positi=(mod(c(3,i+nres),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
+c for each residue check if it is in lipid or lipid water border area
+C respos=mod(c(3,i+nres),boxzsize)
+C print *,positi,bordlipbot,buflipbot
+ if ((positi.gt.bordlipbot)
+ & .and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-
+ &((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i))
+C print *,"I am in true lipid"
+ endif
+ endif ! if in lipid or buffor
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ enddo
+ return
+ end
+
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
SUBROUTINE MATVEC2(A1,V1,V2)
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
return
end
C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+ double precision function sscalelip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+C if(r.lt.r_cut-rlamb) then
+C sscale=1.0d0
+C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C gamm=(r-(r_cut-rlamb))/rlamb
+ sscalelip=1.0d0+r*r*(2*r-3.0d0)
+C else
+C sscale=0d0
+C endif
+ return
+ end
+C-----------------------------------------------------------------------
+ double precision function sscagradlip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+C if(r.lt.r_cut-rlamb) then
+C sscagrad=0.0d0
+C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C gamm=(r-(r_cut-rlamb))/rlamb
+ sscagradlip=r*(6*r-6.0d0)
+C else
+C sscagrad=0.0d0
+C endif
+ return
+ end
+
+C-----------------------------------------------------------------------
double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gvdwpp,
& gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gvdwx,gradcorr,gradxorr,
+ & gliptranc,gliptranx,
& gradcorr5,gradcorr6,gel_loc,gcorr3_turn,gcorr4_turn,gcorr6_turn,
& gel_loc_loc,gel_loc_turn3,gel_loc_turn4,gel_loc_turn6,gcorr_loc,
& g_corr5_loc,g_corr6_loc,gradb,gradbx,gsccorc,gsccorx,gsccor_loc,
& gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres),
& gvdwc(3,maxres),gelc(3,maxres),gvdwpp(3,maxres),
& gradx_scp(3,maxres),
+ & gliptranc(3,-1:maxres),
+ & gliptranx(3,-1:maxres),
& gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres),
& gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres),
& gradcorr5(3,maxres),gradcorr6(3,maxres),
double precision wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
& wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
& wturn6,wvdwpp,wbond,weights,scal14,cutoff_corr,delt_corr,
- & r0_corr
+ & r0_corr,wliptran
integer ipot,n_ene_comp
common /ffield/ wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
& wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
- & wturn6,wvdwpp,wbond,weights(max_ene),
+ & wturn6,wvdwpp,wbond,wliptran,weights(max_ene),
& scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp
common /potentials/ potname(5)
character*3 potname
- double precision aa,bb,augm,aad,bad,app,bpp,ael6,ael3
+ double precision aa_aq,bb_aq,augm,aad,bad,app,bpp,ael6,ael3,
+ & aa_lip,bb_lip
integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro,ielstart,
& ielend,nscp_gr,iscpstart,iscpend,iatsc_s,iatsc_e,iatel_s,
& iatel_e,iatscp_s,iatscp_e,ispp,iscp,expon,expon2
- common /interact/aa(ntyp,ntyp),bb(ntyp,ntyp),augm(ntyp,ntyp),
+ common /interact/aa_aq(ntyp,ntyp),bb_aq(ntyp,ntyp),
+ & augm(ntyp,ntyp),aa_lip(ntyp,ntyp),bb_lip(ntyp,ntyp),
& aad(ntyp,2),bad(ntyp,2),app(2,2),bpp(2,2),ael6(2,2),ael3(2,2),
& expon,expon2,nnt,nct,nint_gr(maxres),istart(maxres,maxint_gr),
& iend(maxres,maxint_gr),itype(maxres),itel(maxres),itypro,
C 12/1/95 Array EPS included in the COMMON block.
double precision eps,sigma,sigmaii,rs0,chi,chip,chip0,alp,signa0,
& sigii,sigma0,rr0,r0,r0e,r0d,rpp,epp,elpp6,elpp3,eps_scp,rscp,
- & eps_orig
+ & eps_orig,epslip
common /body/eps(ntyp,ntyp),sigma(ntyp,ntyp),sigmaii(ntyp,ntyp),
+ &epslip(ntyp,ntyp),
& rs0(ntyp,ntyp),chi(ntyp,ntyp),chip(ntyp),chip0(ntyp),alp(ntyp),
& sigma0(ntyp),sigii(ntyp),rr0(ntyp),r0(ntyp,ntyp),r0e(ntyp,ntyp),
& r0d(ntyp,2),rpp(2,2),epp(2,2),elpp6(2,2),elpp3(2,2),
& aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp),
& distchainmax,nbondterm(ntyp)
&,vbldpDUM
+C 01/29/15 Lipidic parameters
+ double precision pepliptran,liptranene
+ common /lipid/ pepliptran,liptranene(ntyp)
+
ihist=30
iweight=31
izsc=32
+C Lipidic input file for parameters range 60-79
+ iliptranpar=60
C
C Set default weights of the energy terms.
C
enddo
do i=1,ntyp
do j=1,ntyp
- aa(i,j)=0.0D0
- bb(i,j)=0.0D0
+ aa_lip(i,j)=0.0D0
+ bb_lip(i,j)=0.0D0
+ aa_aq(i,j)=0.0D0
+ bb_aq(i,j)=0.0D0
augm(i,j)=0.0D0
sigma(i,j)=0.0D0
r0(i,j)=0.0D0
double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
& eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
- & escloc,
+ & escloc,eliptran,
& ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
& eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
estr=enetb(18,i,iparm)
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
+ eliptran=enetb(22,i,iparm)
#ifdef SPLITELE
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
#else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
#endif
#ifdef MPI
Fdimless(i)=
open (isidep,file=sidename,status='old')
call mygetenv('SIDEP',sidepname)
open (isidep1,file=sidepname,status="old")
+ call mygetenv('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
#ifndef OLDSCP
C
C 8/9/01 In the newest version SCp interaction constants are read from a file
wsccor=ww(19)
whpb=ww(15)
wstrain=ww(15)
+ wliptran=ww(22)
endif
call card_concat(controlcard,.false.)
enddo
enddo
endif
+ read(iliptranpar,*) pepliptran
+ do i=1,ntyp
+ read(iliptranpar,*) liptranene(i)
+ enddo
+ close(iliptranpar)
#ifdef CRYST_THETA
C
C Read the parameters of the probability distribution/energy expression
read (isidep,*)(sigii(i),i=1,ntyp)
read (isidep,*)(chip(i),i=1,ntyp)
read (isidep,*)(alp(i),i=1,ntyp)
+ do i=1,ntyp
+ read (isidep,*)(epslip(i,j),j=i,ntyp)
+C print *,"WARNING!!"
+C do j=1,ntyp
+C epslip(i,j)=epslip(i,j)+0.05d0
+C enddo
+ enddo
C For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do i=1,ntyp
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
do i=1,ntyp
do j=i,ntyp
epsij=eps(i,j)
+ epsijlip=epslip(i,j)
if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
rrij=sigma(i,j)
else
epsij=eps(i,j)
sigeps=dsign(1.0D0,epsij)
epsij=dabs(epsij)
- aa(i,j)=epsij*rrij*rrij
- bb(i,j)=-sigeps*epsij*rrij
- aa(j,i)=aa(i,j)
- bb(j,i)=bb(i,j)
+ aa_aq(i,j)=epsij*rrij*rrij
+ bb_aq(i,j)=-sigeps*epsij*rrij
+ aa_aq(j,i)=aa_aq(i,j)
+ bb_aq(j,i)=bb_aq(i,j)
+ sigeps=dsign(1.0D0,epsijlip)
+ epsijlip=dabs(epsijlip)
+ aa_lip(i,j)=epsijlip*rrij*rrij
+ bb_lip(i,j)=-sigeps*epsijlip*rrij
+ aa_lip(j,i)=aa_lip(i,j)
+ bb_lip(j,i)=bb_lip(i,j)
if (ipot.gt.2) then
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
c Cutoff range for interactions
call reada(controlcard,"R_CUT",r_cut,15.0d0)
call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+ call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+ call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+ if (lipthick.gt.0.0d0) then
+ bordliptop=(boxzsize+lipthick)/2.0
+ bordlipbot=bordliptop-lipthick
+C endif
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
+ & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+ buflipbot=bordlipbot+lipbufthick
+ bufliptop=bordliptop-lipbufthick
+ if ((lipbufthick*2.0d0).gt.lipthick)
+ &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+ endif
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
call readi(controlcard,'SYM',symetr,1)
write (iout,*) "DISTCHAINMAX",distchainmax
write (iout,*) "delta",delta
ww_all(16,iparm)=wvdwpp
ww_all(17,iparm)=wbond
ww_all(19,iparm)=wsccor
+ ww_all(22,iparm)=wliptran
c Store bond parameters
vbldp0_all(iparm)=vbldp0
akp_all(iparm)=akp
c Store sidechain parameters
do i=1,ntyp
do j=1,ntyp
- aa_all(j,i,iparm)=aa(j,i)
- bb_all(j,i,iparm)=bb(j,i)
+ aa_aq_all(j,i,iparm)=aa_aq(j,i)
+ bb_aq_all(j,i,iparm)=bb_aq(j,i)
+ aa_lip_all(j,i,iparm)=aa_lip(j,i)
+ bb_lip_all(j,i,iparm)=bb_lip(j,i)
r0_all(j,i,iparm)=r0(j,i)
sigma_all(j,i,iparm)=sigma(j,i)
chi_all(j,i,iparm)=chi(j,i)
augm_all(j,i,iparm)=augm(j,i)
eps_all(j,i,iparm)=eps(j,i)
+ epslip_all(j,i,iparm)=epslip(j,i)
enddo
enddo
do i=1,ntyp
wvdwpp=ww_all(16,iparm)
wbond=ww_all(17,iparm)
wsccor=ww_all(19,iparm)
+ wliptran=ww_all(22,iparm)
c Restore bond parameters
vbldp0=vbldp0_all(iparm)
akp=akp_all(iparm)
c Restore sidechain parameters
do i=1,ntyp
do j=1,ntyp
- aa(j,i)=aa_all(j,i,iparm)
- bb(j,i)=bb_all(j,i,iparm)
+ aa_aq(j,i)=aa_aq_all(j,i,iparm)
+ bb_aq(j,i)=bb_aq_all(j,i,iparm)
+ aa_lip(j,i)=aa_lip_all(j,i,iparm)
+ bb_lip(j,i)=bb_lip_all(j,i,iparm)
r0(j,i)=r0_all(j,i,iparm)
sigma(j,i)=sigma_all(j,i,iparm)
chi(j,i)=chi_all(j,i,iparm)
augm(j,i)=augm_all(j,i,iparm)
eps(j,i)=eps_all(j,i,iparm)
+ epslip(j,i)=epslip_all(j,i,iparm)
enddo
enddo
do i=1,ntyp
& eplus,eminus,logfac,tanhT,tt
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
& escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
- & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+ & eliptran
integer ind_point(maxpoint),upindE,indE
character*16 plik
estr=enetb(18,i,iparm)
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
+ eliptran=enetb(22,i,iparm)
+
#ifdef DEBUG
write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
& evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
#else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
#endif
#ifdef DEBUG
write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
& +ftprim(1)*wtor*etors+
& ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
& +ftprim(1)*wtor*etors+
& ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
projects / unres.git / commitdiff