X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=9271240aaff989920105157fb8a10df0578752c8;hb=d7d33ac73e7386cf74bb07da09b67fa6619a0b11;hp=346a3c3200b1acae0d76e979343576bac8e5bb39;hpb=befb4994bdeb1dd3bb9d17881522e33b31bc0ea1;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 346a3c3..9271240 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -99,6 +99,7 @@ c endif C C Compute the side-chain and electrostatic interaction energy C +C print *,ipot goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj(evdw) @@ -112,6 +113,7 @@ C Berne-Pechukas potential (dilated LJ, angular dependence). goto 107 C Gay-Berne potential (shifted LJ, angular dependence). 104 call egb(evdw) +C print *,"bylem w egb" goto 107 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). 105 call egbv(evdw) @@ -122,6 +124,11 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +cmc +cmc Sep-06: egb takes care of dynamic ss bonds too +cmc +c if (dyn_ss) call dyn_set_nss + c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time01=MPI_Wtime() @@ -152,9 +159,9 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" eello_turn4=0.0d0 endif else -c write (iout,*) "Soft-spheer ELEC potential" - call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, - & eello_turn4) + write (iout,*) "Soft-spheer ELEC potential" +c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, +c & eello_turn4) endif c print *,"Processor",myrank," computed UELEC" C @@ -194,6 +201,7 @@ c print *,"Processor",myrank," computed UB" C C Calculate the SC local energy. C +C print *,"TU DOCHODZE?" call esc(escloc) c print *,"Processor",myrank," computed USC" C @@ -224,6 +232,7 @@ C else esccor=0.0d0 endif +C print *,"PRZED MULIt" c print *,"Processor",myrank," computed Usccorr" C C 12/1/95 Multi-body terms @@ -256,6 +265,14 @@ C after the equilibration time Uconst=0.0d0 Uconst_back=0.0d0 endif +C 01/27/2015 added by adasko +C the energy component below is energy transfer into lipid environment +C based on partition function +C print *,"przed lipidami" + if (wliptran.gt.0) then + call Eliptransfer(eliptran) + endif +C print *,"za lipidami" #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -297,10 +314,12 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor + energia(22)=eliptran c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) + if (dyn_ss) call dyn_set_nss c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 @@ -387,20 +406,21 @@ cMS$ATTRIBUTES C :: proc_proc 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 c detecting NaNQ @@ -436,9 +456,10 @@ cMS$ATTRIBUTES C :: proc_proc #endif #ifdef MPI include 'mpif.h' - double precision gradbufc(3,maxres),gradbufx(3,maxres), - & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) #endif + double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres), + & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres) + & ,gloc_scbuf(3,-1:maxres) include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.FFIELD' @@ -491,7 +512,7 @@ c enddo 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))+ @@ -502,10 +523,12 @@ c enddo & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) + & +wliptran*gliptranc(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))+ @@ -517,6 +540,7 @@ c enddo & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) + & +wliptran*gliptranc(j,i) enddo enddo #endif @@ -530,7 +554,7 @@ c enddo 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 @@ -573,7 +597,7 @@ c enddo 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 @@ -594,7 +618,7 @@ c enddo enddo call flush(iout) #endif - do i=1,nres + do i=-1,nres do j=1,3 gradbufc_sum(j,i)=gradbufc(j,i) gradbufc(j,i)=0.0d0 @@ -603,7 +627,7 @@ c enddo 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 @@ -631,7 +655,7 @@ c enddo do k=1,3 gradbufc(k,nres)=0.0d0 enddo - 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)+ @@ -652,6 +676,7 @@ c enddo & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) + & +wliptran*gliptranc(j,i) #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ @@ -671,12 +696,14 @@ c enddo & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) + & +wliptran*gliptranc(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) + & +wliptran*gliptranx(j,i) enddo enddo #ifdef DEBUG @@ -965,6 +992,7 @@ C------------------------------------------------------------------------ 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, @@ -973,7 +1001,7 @@ C------------------------------------------------------------------------ & 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)'/ @@ -997,6 +1025,7 @@ C------------------------------------------------------------------------ & '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, @@ -1005,7 +1034,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,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)'/ @@ -1028,6 +1057,7 @@ C------------------------------------------------------------------------ & '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 @@ -1082,13 +1112,14 @@ C Change 12/1/95 to calculate four-body interactions c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C have you changed here? + e1=fac*fac*aa + e2=fac*bb evdwij=e1+e2 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') -cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), +cd & restyp(itypi),i,restyp(itypj),j,a(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) evdw=evdw+evdwij @@ -1232,8 +1263,9 @@ C rij=1.0D0/r_inv_ij r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C have you changed here? + e1=fac*fac*aa + e2=fac*bb evdwij=e_augm+e1+e2 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -1359,17 +1391,18 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. call sc_angular C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives +C have you changed here? 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 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/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa cd write (iout,'(2(a3,i3,2x),15(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & epsi,sigm,chi1,chi2,chip1,chip2, @@ -1414,20 +1447,21 @@ C include 'COMMON.CALC' include 'COMMON.CONTROL' include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' logical lprn integer xshift,yshift,zshift evdw=0.0D0 ccccc energy_dec=.false. -c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon +C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0 C we have the original box) - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle @@ -1466,6 +1500,38 @@ c endif if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize +C define scaling factor for lipids + +C 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 ((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 + +C xi=xi+xshift*boxxsize +C yi=yi+yshift*boxysize +C zi=zi+zshift*boxzsize dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) @@ -1479,6 +1545,12 @@ C Calculate SC interaction energy. C 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' + ELSE ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1541,12 +1613,73 @@ c endif 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 (zi.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 + if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)') + &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) +C if (ssgradlipj.gt.0.0d0) print *,"??WTF??" + print *,sslipi,sslipj,bordlipbot,zi,zj + 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 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) - xj=xj-xi - yj=yj-yi - zj=zj-zi +C xj=xj-xi +C yj=yj-yi +C zj=zj-zi c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi c write (iout,*) "j",j," dc_norm", c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) @@ -1578,18 +1711,24 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C here to start with +C if (c(i,3).gt. + faclip=fac + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt +C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij, +C &((sslipi+sslipj)/2.0d0+ +C &(2.0d0-sslipi-sslipj)/2.0d0) c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 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/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@ -1611,18 +1750,27 @@ c & evdwij,fac,sigma(itypi,itypj),expon fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij c fac=0.0d0 C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt) + &*(eps3rt*eps3rt)*sss/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 +C gg_lipi(3)=0.0d0 +C gg_lipj(3)=0.0d0 gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac C Calculate angular part of the gradient. call sc_grad endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i - enddo ! zshift - enddo ! yshift - enddo ! xshift +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return @@ -1659,6 +1807,41 @@ c if (icall.eq.0) lprn=.true. xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + xi=mod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=mod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize +C define scaling factor for lipids + +C 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 ((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- + & ((positi-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-positi)/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) @@ -1695,9 +1878,74 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi +C xj=c(1,nres+j)-xi +C yj=c(2,nres+j)-yi +C zj=c(3,nres+j)-zi + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + 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- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/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 if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5') +C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) + 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 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1718,8 +1966,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -1728,8 +1976,8 @@ c--------------------------------------------------------------- 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/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa 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), @@ -1743,6 +1991,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac-2*expon*rrij*e_augm + fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -1853,10 +2102,10 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 enddo c 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 - 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 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) @@ -1873,8 +2122,8 @@ cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) cgrad enddo cgrad 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 @@ -1976,7 +2225,7 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' dimension ggg(3) -cd write(iout,*) 'In EELEC_soft_sphere' +C write(iout,*) 'In EELEC_soft_sphere' ees=0.0D0 evdw1=0.0D0 eel_loc=0.0d0 @@ -1991,6 +2240,12 @@ cd write(iout,*) 'In EELEC_soft_sphere' xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -2004,10 +2259,49 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) dxj=dc(1,j) dyj=dc(2,j) dzj=dc(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj + xj=mod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=mod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif rij=xj*xj+yj*yj+zj*zj + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) if (rij.lt.r0ijsq) then evdw1ij=0.25d0*(rij-r0ijsq)**2 fac=rij-r0ijsq @@ -2015,13 +2309,13 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) evdw1ij=0.0d0 fac=0.0d0 endif - evdw1=evdw1+evdw1ij + evdw1=evdw1+evdw1ij*sss C C Calculate contributions to the Cartesian gradient. C - ggg(1)=fac*xj - ggg(2)=fac*yj - ggg(3)=fac*zj + ggg(1)=fac*xj*sssgrad + ggg(2)=fac*yj*sssgrad + ggg(3)=fac*zj*sssgrad do k=1,3 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) @@ -2879,6 +3173,8 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups C C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end + if (i.le.1) cycle +C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 & .or. itype(i+3).eq.ntyp1 @@ -2894,31 +3190,6 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi -C Return atom into box, boxxsize is size of box in x dimension -c 184 continue -c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize -c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize -C Condition for being inside the proper box -c if ((xmedi.gt.((0.5d0)*boxxsize)).or. -c & (xmedi.lt.((-0.5d0)*boxxsize))) then -c go to 184 -c endif -c 185 continue -c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize -c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize -cC Condition for being inside the proper box -c if ((ymedi.gt.((0.5d0)*boxysize)).or. -c & (ymedi.lt.((-0.5d0)*boxysize))) then -c go to 185 -c endif -c 186 continue -c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize -c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize -cC Condition for being inside the proper box -c if ((zmedi.gt.((0.5d0)*boxzsize)).or. -c & (zmedi.lt.((-0.5d0)*boxzsize))) then -c go to 186 -c endif xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) @@ -2931,6 +3202,7 @@ c endif num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end + if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 @@ -2986,13 +3258,14 @@ c endif num_cont_hb(i)=num_conti enddo ! i C Loop over all neighbouring boxes - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e + if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 & .or. itype(i-1).eq.ntyp1 @@ -3012,8 +3285,11 @@ c if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize +C xmedi=xmedi+xshift*boxxsize +C ymedi=ymedi+yshift*boxysize +C zmedi=zmedi+zshift*boxzsize -C Return atom into box, boxxsize is size of box in x dimension +C Return tom into box, boxxsize is size of box in x dimension c 164 continue c if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize c if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize @@ -3042,7 +3318,8 @@ c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) -c write (iout,*) i,j,itype(i),itype(j) +C write (iout,*) i,j + if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 & .or.itype(j+2).eq.ntyp1 & .or.itype(j-1).eq.ntyp1 @@ -3051,9 +3328,9 @@ c write (iout,*) i,j,itype(i),itype(j) enddo ! j num_cont_hb(i)=num_conti enddo ! i - enddo ! zshift - enddo ! yshift - enddo ! xshift +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift c write (iout,*) "Number of loop steps in EELEC:",ind cd do i=1,nres @@ -3132,7 +3409,38 @@ C zj=c(3,j)+0.5D0*dzj-zmedi if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize - + if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ" + dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + isubchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + isubchap=1 + endif + enddo + enddo + enddo + if (isubchap.eq.1) then + xj=xj_temp-xmedi + yj=yj_temp-ymedi + zj=zj_temp-zmedi + else + xj=xj_safe-xmedi + yj=yj_safe-ymedi + zj=zj_safe-zmedi + endif C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC c 174 continue c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize @@ -3159,9 +3467,9 @@ c & (zj.lt.((-0.5d0)*boxzsize))) then c go to 176 c endif C endif !endPBC condintion - xj=xj-xmedi - yj=yj-ymedi - zj=zj-zmedi +C xj=xj-xmedi +C yj=yj-ymedi +C zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj sss=sscale(sqrt(rij)) @@ -4104,9 +4412,9 @@ C r0_scp=4.5d0 cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) @@ -4144,7 +4452,9 @@ c endif if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - +C xi=xi+xshift*boxxsize +C yi=yi+yshift*boxysize +C zi=zi+zshift*boxzsize do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -4187,10 +4497,41 @@ c go to 176 if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize + 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 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif c c endif - xj=xj-xi - yj=yj-yi - zj=zj-zi +C xj=xj-xi +C yj=yj-yi +C zj=zj-zi rij=xj*xj+yj*yj+zj*zj r0ij=r0_scp @@ -4243,9 +4584,9 @@ cgrad enddo enddo ! iint enddo ! i - enddo !zshift - enddo !yshift - enddo !xshift +C enddo !zshift +C enddo !yshift +C enddo !xshift return end C----------------------------------------------------------------------------- @@ -4270,11 +4611,12 @@ C dimension ggg(3) evdw2=0.0D0 evdw2_14=0.0d0 +c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) @@ -4287,7 +4629,10 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize - +c xi=xi+xshift*boxxsize +c yi=yi+yshift*boxysize +c zi=zi+zshift*boxzsize +c print *,xi,yi,zi,'polozenie i' C Return atom into box, boxxsize is size of box in x dimension c 134 continue c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize @@ -4298,6 +4643,8 @@ c & (xi.lt.((xshift-0.5d0)*boxxsize))) then c go to 134 c endif c 135 continue +c print *,xi,boxxsize,"pierwszy" + c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize C Condition for being inside the proper box @@ -4356,13 +4703,46 @@ c if ((zj.gt.((0.5d0)*boxzsize)).or. c & (zj.lt.((-0.5d0)*boxzsize))) then c go to 176 c endif - xj=xj-xi - yj=yj-yi - zj=zj-zi +CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE + 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 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif +c print *,xj,yj,zj,'polozenie j' rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +c print *,rrij sss=sscale(1.0d0/(dsqrt(rrij))) +c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz' +c if (sss.eq.0) print *,'czasem jest OK' + if (sss.le.0.0d0) cycle sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) - if (sss.gt.0.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) @@ -4415,14 +4795,14 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - endif !endif for sscale cutoff +c endif !endif for sscale cutoff enddo ! j enddo ! iint enddo ! i - enddo !zshift - enddo !yshift - enddo !xshift +c enddo !zshift +c enddo !yshift +c enddo !xshift do i=1,nct do j=1,3 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) @@ -4472,50 +4852,58 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif -cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj +c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, +c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. - & iabs(itype(jjj)).eq.1) then +C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C & iabs(itype(jjj)).eq.1) then +cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds + if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond - additional check + if (ii.gt.nres + & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij else C Calculate the distance between the two points and its difference from the C target distance. - dd=dist(ii,jj) - rdis=dd-dhpb(i) + dd=dist(ii,jj) + rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. - waga=forcon(i) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis C C Evaluate gradient. C - fac=waga*rdis/dd + fac=waga*rdis/dd cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). - if (iii.lt.ii) then + if (iii.lt.ii) then do j=1,3 ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo - endif + endif cgrad do j=iii,jjj-1 cgrad do k=1,3 cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k) cgrad enddo cgrad enddo - do k=1,3 - ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) - ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) - enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo endif enddo ehpb=0.5D0*ehpb @@ -4648,7 +5036,7 @@ C YES vbldpDUM is the equlibrium length of spring for Dummy atom C NO vbldp0 is the equlibrium lenght of spring for peptide group diff = vbld(i)-vbldp0 endif - if (energy_dec) write (iout,'(a7,i5,4f7.3)') + if (energy_dec) write (iout,'(a7,i5,4f7.3)') & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 @@ -4667,7 +5055,7 @@ c nbi=nbondterm(iti) if (nbi.eq.1) then diff=vbld(i+nres)-vbldsc0(1,iti) - if (energy_dec) write (iout,*) + if (energy_dec) write (iout,*) & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, & AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff @@ -8475,12 +8863,12 @@ C C C o o C C \ /l\ /j\ / C C \ / \ / \ / C -C o| o | | o |o C +C o| o | | o |o C C \ j|/k\| \ |/k\|l C C \ / \ \ / \ C C o o C -C i i C -C C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l C AL 7/4/01 s1 would occur in the sixth-order moment, @@ -8651,10 +9039,10 @@ c---------------------------------------------------------------------------- double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C C +C C C Parallel Antiparallel C C C -C o o C +C o o C C /l\ / \ /j\ C C / \ / \ / \ C C /| o |o o| o |\ C @@ -8768,7 +9156,7 @@ c---------------------------------------------------------------------------- & auxvec1(2),auxmat1(2,2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C C +C C C Parallel Antiparallel C C C C o o C @@ -8776,10 +9164,10 @@ C /l\ / \ /j\ C C / \ / \ / \ C C /| o |o o| o |\ C C \ j|/k\| \ |/k\|l C -C \ / \ \ / \ C +C \ / \ \ / \ C C o \ o \ C C i i C -C C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective @@ -9479,4 +9867,120 @@ crc print *,((prod(i,j),i=1,2),j=1,2) return end +CCC---------------------------------------------- + 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.NAMES' + 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=ilip_start,ilip_end +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 + +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*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" + 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 +C do i=1,1 + 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))/2.0d0 + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i))/2.0d0 +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))/2.0d0 + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i))/2.0d0 +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