X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=7207b35c4f2b5075be6cefe490c0accd9508166f;hp=a8c0323eb73e7ad588b76621e6ba7f00c3e9b221;hb=0c5fb524ca994d617a32adfb678302e2f3dcef1e;hpb=ffd0e6dcc64b1ef30773e02b4ef72c1b2455fea5 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 a8c0323..7207b35 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -199,6 +199,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 @@ -229,6 +230,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 @@ -264,9 +266,11 @@ C after the equilibration time 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 + call Eliptransfer(eliptran) endif +C print *,"za lipidami" #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -308,7 +312,7 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor - energia(22)=eliptrans + 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" @@ -400,7 +404,7 @@ cMS$ATTRIBUTES C :: proc_proc estr=energia(17) Uconst=energia(20) esccor=energia(21) - energia(22)=eliptrans + eliptran=energia(22) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc @@ -451,8 +455,8 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI include 'mpif.h' #endif - double precision gradbufc(3,maxres),gradbufx(3,maxres), - & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) + double precision gradbufc(3,0:maxres),gradbufx(3,0:maxres), + & glocbuf(4*maxres),gradbufc_sum(3,0:maxres),gloc_scbuf(3,0:maxres) include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.FFIELD' @@ -516,6 +520,8 @@ c enddo & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) + & +wliptran*gliptranc(j,i) + enddo enddo #else @@ -531,6 +537,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 @@ -982,6 +989,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, @@ -990,7 +998,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)'/ @@ -1014,6 +1022,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, @@ -1022,7 +1031,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)'/ @@ -1045,6 +1054,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 @@ -1099,6 +1109,7 @@ 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 +C have you changed here? e1=fac*fac*aa(itypi,itypj) e2=fac*bb(itypi,itypj) evdwij=e1+e2 @@ -1249,6 +1260,7 @@ C rij=1.0D0/r_inv_ij r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon +C have you changed here? e1=fac*fac*aa(itypi,itypj) e2=fac*bb(itypi,itypj) evdwij=e_augm+e1+e2 @@ -1376,6 +1388,7 @@ 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) @@ -1636,6 +1649,8 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon +C here to start with +C if (c(i,3).gt. e1=fac*fac*aa(itypi,itypj) e2=fac*bb(itypi,itypj) evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1674,6 +1689,7 @@ C Calculate the radial part of the gradient gg(3)=zj*fac C Calculate angular part of the gradient. call sc_grad + endif ENDIF ! dyn_ss enddo ! j enddo ! iint @@ -2982,6 +2998,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 @@ -3009,6 +3027,7 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition 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 @@ -3071,6 +3090,7 @@ 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 @@ -3123,7 +3143,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 @@ -4660,15 +4681,14 @@ 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 ->>>>>>> f5379d3246c4bd95e946c4d35d4a1c13e329c4cb call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij endif @@ -9674,6 +9694,7 @@ crc print *,((prod(i,j),i=1,2),j=1,2) end CCC---------------------------------------------- subroutine Eliptransfer(eliptran) + implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -9687,6 +9708,7 @@ CCC---------------------------------------------- include 'COMMON.CONTROL' include 'COMMON.SPLITELE' include 'COMMON.SBRIDGE' +C print *,"wchodze" C structure of box: C water C--bordliptop-- buffore starts @@ -9695,67 +9717,91 @@ C lipid C--buflipbot--- lipid ends buffore starts C--bordlipbot--buffore ends eliptran=0.0 - do i=1,nres + do i=ilip_start,ilip_end + if (itype(i).eq.ntyp1) cycle + + positi=(mod((c(3,i)+c(3,i+1)),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 ((mod(c(3,i),boxzsize).gt.bordlipbot) - &.and.(mod(c(3,i),boxzsize).lt.bordliptop)) then + if ((positi.gt.bordlipbot) + &.and.(positi.lt.bordliptop)) then C the energy transfer exist - if (mod(c(3,i),boxzsize).lt.buflipbot) then + if (positi.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- - & ((mod(c(3,i),boxzsize)-bordlipbot)/lipbufthick) + & ((positi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore - ssslip=sscale(fracinbuf) + sslip=sscalelip(fracinbuf) ssgradlip=-sscagradlip(fracinbuf)/lipbufthick - eliptran=eliptran+sslip - gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0 +C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran + C print *,"doing sccale for lower part" - elseif (mod(c(3,i),boxzsize).gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-mod(c(3,i),boxzsize))/lipbufthick) - ssslip=sscale(fracinbuf) + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) ssgradlip=sscagradlip(fracinbuf)/lipbufthick - eliptran=eliptran+sslip - gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran + 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 print *, "doing sscalefor top part" else - eliptran=eliptran+1.0d0 + eliptran=eliptran+pepliptran 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 - eliptran=eliptran*pepliptran +C eliptran=eliptran*pepliptran C now the same for side chains - do i=1,nres + 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 - if ((mod(c(3,i+nres),boxzsize).gt.bordlipbot) - & .and.(mod(c(3,i+nres),boxzsize).lt.bordliptop)) then +C respos=mod(c(3,i+nres),boxzsize) + if ((positi.gt.bordlipbot) + & .and.(positi.lt.bordliptop)) then C the energy transfer exist - if (mod(c(3,i+nres),boxzsize).lt.buflipbot) then + if (positi.lt.buflipbot) then fracinbuf=1.0d0- - & ((mod(c(3,i+nres),boxzsize)-bordlipbot)/lipbufthick) + & ((positi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore - ssslip=sscale(fracinbuf) + sslip=sscalelip(fracinbuf) ssgradlip=-sscagradlip(fracinbuf)/lipbufthick eliptran=eliptran+sslip*liptranene(itype(i)) - gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i))/2.0d0 + gliptranc(3,i-1)= + &+ssgradlip*liptranene(itype(i)) print *,"doing sccale for lower part" - elseif (mod(c(3,i+nres),boxzsize).gt.bufliptop) then + elseif (positi.gt.bufliptop) then fracinbuf=1.0d0- - &((bordliptop-mod(c(3,i+nres),boxzsize))/lipbufthick) - ssslip=sscale(fracinbuf) + &((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)) - print *, "doing sscalefor top part" + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i))/2.0d0 + gliptranc(3,i-1)= + &+ssgradlip*liptranene(itype(i)) + print *, "doing sscalefor top part",sslip,fracinbuf else eliptran=eliptran+liptranene(itype(i)) 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