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=e6f630d1303b8d4ab6a7acd72c817eda74d05333;hb=29fc8e85b7adc50ed8c61c03036da4e08cddfedd;hp=c13d3410c2a4813d243821240e518f09072d0855;hpb=839e87c798b8f0463a18e3ea4bf5105108546fed;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 c13d341..e6f630d 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -24,6 +24,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.MD' include 'COMMON.CONTROL' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -98,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) @@ -111,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) @@ -121,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() @@ -151,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 @@ -193,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 @@ -223,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 @@ -255,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 @@ -296,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 @@ -386,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 @@ -435,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' @@ -490,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))+ @@ -501,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))+ @@ -516,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 @@ -529,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 @@ -572,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 @@ -593,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 @@ -602,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 @@ -630,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)+ @@ -651,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)+ @@ -670,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 @@ -710,7 +738,7 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -719,7 +747,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG do i=1,nres do j=1,3 gloc_scbuf(j,i)=gloc_sc(j,i,icg) @@ -739,7 +767,7 @@ c enddo call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -748,7 +776,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -964,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, @@ -972,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)'/ @@ -996,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, @@ -1004,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)'/ @@ -1027,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 @@ -1081,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 @@ -1231,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) @@ -1358,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, @@ -1412,14 +1446,22 @@ C include 'COMMON.IOUNITS' 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) +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 @@ -1427,6 +1469,70 @@ c if (icall.eq.0) lprn=.false. xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+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 +c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize +C Condition for being inside the proper box +c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. +c & (xi.lt.((xshift-0.5d0)*boxxsize))) then +c go to 134 +c endif +c 135 continue +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 +c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. +c & (yi.lt.((yshift-0.5d0)*boxysize))) then +c go to 135 +c endif +c 136 continue +c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize +c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize +C Condition for being inside the proper box +c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 136 +c endif + 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- + & ((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) dzi=dc_norm(3,nres+i) @@ -1439,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 @@ -1467,17 +1579,118 @@ 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 + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) +C Return atom J into box the original box +c 137 continue +c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize +c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box +c if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 137 +c endif +c 138 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +C Condition for being inside the proper box +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 138 +c endif +c 139 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +C Condition for being inside the proper box +c if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 139 +c endif + 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- + & ((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 +C if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)') +C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) +C if (ssgradlipj.gt.0.0d0) print *,"??WTF??" +C 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) +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) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) + sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) + +c write (iout,'(a7,4f8.3)') +c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb + if (sss.gt.0.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -1498,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 + 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, @@ -1526,16 +1745,32 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac +c print '(2i4,6f8.4)',i,j,sss,sssgrad* +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 +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return @@ -1572,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) @@ -1608,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) @@ -1631,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 @@ -1641,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), @@ -1656,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 @@ -1743,6 +2079,7 @@ C---------------------------------------------------------------------------- include 'COMMON.CALC' include 'COMMON.IOUNITS' double precision dcosom1(3),dcosom2(3) +cc print *,'sss=',sss eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 @@ -1761,16 +2098,16 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k)) enddo do k=1,3 - gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k) + gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss 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 - gvdwx(k,j)=gvdwx(k,j)+gg(k) + & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss + 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 + & +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)) c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) @@ -1785,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 @@ -1888,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 @@ -1903,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) @@ -1916,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 @@ -1927,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) @@ -2252,7 +2634,7 @@ C C Compute the virtual-bond-torsional-angle dependent quantities needed C to calculate the el-loc multibody terms of various order. C - write(iout,*) 'nphi=',nphi,nres +c write(iout,*) 'nphi=',nphi,nres #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else @@ -2270,23 +2652,23 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then else iti1=ntortyp+1 endif - write(iout,*),i - b1(1,i-2)=bnew1(1,1,iti)*sin(theta(i-1)/2.0) - & +bnew1(2,1,iti)*sin(theta(i-1)) - & +bnew1(3,1,iti)*cos(theta(i-1)/2.0) - gtb1(1,i-2)=bnew1(1,1,iti)*cos(theta(i-1)/2.0)/2.0 - & +bnew1(2,1,iti)*cos(theta(i-1)) - & -bnew1(3,1,iti)*sin(theta(i-1)/2.0)/2.0 +c write(iout,*),i + b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0) + & +bnew1(2,1,iti)*dsin(theta(i-1)) + & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0) + gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 + & +bnew1(2,1,iti)*dcos(theta(i-1)) + & -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) - b2(1,i-2)=bnew2(1,1,iti)*sin(theta(i-1)/2.0) - & +bnew2(2,1,iti)*sin(theta(i-1)) - & +bnew2(3,1,iti)*cos(theta(i-1)/2.0) + b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0) + & +bnew2(2,1,iti)*dsin(theta(i-1)) + & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0) c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) - gtb2(1,i-2)=bnew2(1,1,iti)*cos(theta(i-1)/2.0)/2.0 - & +bnew2(2,1,iti)*cos(theta(i-1)) - & -bnew2(3,1,iti)*sin(theta(i-1)/2.0)/2.0 + gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 + & +bnew2(2,1,iti)*dcos(theta(i-1)) + & -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 c if (ggb1(1,i).eq.0.0d0) then c write(iout,*) 'i=',i,ggb1(1,i), c &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0, @@ -2297,7 +2679,14 @@ c endif gtb1(2,i-2)=0.0 b2(2,i-2)=bnew2(1,2,iti) gtb2(2,i-2)=0.0 -c EE(1,1,iti)=0.0d0 + EE(1,1,i-2)=eenew(1,iti)*dcos(theta(i-1)) + EE(1,2,i-2)=eeold(1,2,iti) + EE(2,1,i-2)=eeold(2,1,iti) + EE(2,2,i-2)=eeold(2,2,iti) + gtEE(1,1,i-2)=-eenew(1,iti)*dsin(theta(i-1)) + gtEE(1,2,i-2)=0.0d0 + gtEE(2,2,i-2)=0.0d0 + gtEE(2,1,i-2)=0.0d0 c EE(2,2,iti)=0.0d0 c EE(1,2,iti)=0.5d0*eenew(1,iti) c EE(2,1,iti)=0.5d0*eenew(1,iti) @@ -2307,8 +2696,9 @@ c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) b1tilde(2,i-2)=-b1(2,i-2) b2tilde(1,i-2)=b2(1,i-2) b2tilde(2,i-2)=-b2(2,i-2) - write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) - write (iout,*) 'theta=', theta(i-1) +c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) +c write(iout,*) 'b1=',b1(1,i-2) +c write (iout,*) 'theta=', theta(i-1) enddo #ifdef PARMAT do i=ivec_start+2,ivec_end+2 @@ -2383,19 +2773,17 @@ c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) Ug2der(2,2,i-2)=0.0d0 endif c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then -#ifndef NEWCORR if (i.gt. nnt+2 .and. i.lt.nct+2) then iti = itortyp(itype(i-2)) else - iti=ntortyp+1 + iti=ntortyp 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 + iti1=ntortyp endif -#endif cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) cd write (iout,*) 'b2',b2(:,iti) @@ -2407,7 +2795,13 @@ c if (i .gt. iatel_s+2) then call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2)) c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" #endif - call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2)) +c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti), +c & EE(1,2,iti),EE(2,2,iti) + call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2)) + call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2)) +c write(iout,*) "Macierz EUG", +c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2), +c & eug(2,2,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2)) @@ -2430,7 +2824,7 @@ c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" enddo endif call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2)) - call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2)) + call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2)) do k=1,2 muder(k,i-2)=Ub2der(k,i-2) enddo @@ -2439,15 +2833,15 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (itype(i-1).le.ntyp) then iti1 = itortyp(itype(i-1)) else - iti1=ntortyp+1 + iti1=ntortyp endif else - iti1=ntortyp+1 + iti1=ntortyp endif do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo -cd write (iout,*) 'mu ',mu(:,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) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) @@ -2771,6 +3165,7 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -2853,9 +3248,16 @@ c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms C 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) cycle + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2865,16 +3267,26 @@ C 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 TU ZLE -c call eelecij(i,i+2,ees,evdw1,eel_loc) + 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 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) cycle + & .or. itype(i+4).eq.ntyp1 + & .or. itype(i+5).eq.ntyp1 + & .or. itype(i).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2884,19 +3296,58 @@ c call eelecij(i,i+2,ees,evdw1,eel_loc) 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 194 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 194 +c endif +c 195 continue +c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize +c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize +C 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 195 +c endif +c 196 continue +c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize +c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +C 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 196 +c endif + 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=num_cont_hb(i) -c TU ZLE -c call eelecij(i,i+3,ees,evdw1,eel_loc) +c write(iout,*) "JESTEM W PETLI" + call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti enddo ! i +C Loop over all neighbouring boxes +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 -c do i=iatel_s,iatel_e - do i=4,5 - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + 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 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2906,16 +3357,59 @@ c do i=iatel_s,iatel_e 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 +C xmedi=xmedi+xshift*boxxsize +C ymedi=ymedi+yshift*boxysize +C zmedi=zmedi+zshift*boxzsize + +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 +C Condition for being inside the proper box +c if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or. +c & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then +c go to 164 +c endif +c 165 continue +c if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize +c if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize +C Condition for being inside the proper box +c if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or. +c & (ymedi.lt.((yshift-0.5d0)*boxysize))) then +c go to 165 +c endif +c 166 continue +c if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize +c if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +cC Condition for being inside the proper box +c if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 166 +c endif + c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) -c do j=ielstart(i),ielend(i) - do j=8,9 - write (iout,*) 'tu wchodze',i,j,itype(i),itype(j) - if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle + do j=ielstart(i),ielend(i) +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 + &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti enddo ! i +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift + c write (iout,*) "Number of loop steps in EELEC:",ind cd do i=1,nres cd write (iout,'(i3,3f10.5,5x,3f10.5)') @@ -2946,6 +3440,7 @@ C------------------------------------------------------------------------------- include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -2981,10 +3476,84 @@ c ind=ind+1 dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(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 +C xj=c(1,j)+0.5D0*dxj-xmedi +C yj=c(2,j)+0.5D0*dyj-ymedi +C 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 + 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 +c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box +c if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 174 +c endif +c 175 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +C Condition for being inside the proper box +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 175 +c endif +c 176 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +C Condition for being inside the proper box +c if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 176 +c endif +C endif !endPBC condintion +C xj=xj-xmedi +C yj=yj-ymedi +C zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj + + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) +c if (sss.gt.0.0d0) then rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij @@ -3000,14 +3569,15 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions ev2=bbb*r6ij fac3=ael6i*r6ij fac4=ael3i*r3ij - evdwij=ev1+ev2 + evdwij=(ev1+ev2) el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac - eesij=el1+el2 +C MARYSIA + eesij=(el1+el2) C 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 - evdw1=evdw1+evdwij + evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, @@ -3024,7 +3594,7 @@ C C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij) + facvdw=-6*rrmij*(ev1+evdwij)*sss facel=-3*rrmij*(el1+eesij) fac1=fac erij(1)=xj*rmij @@ -3054,9 +3624,15 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + if (sss.gt.0.0) then + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + else + ggg(1)=0.0 + ggg(2)=0.0 + ggg(3)=0.0 + endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -3076,8 +3652,9 @@ cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) cgrad enddo cgrad enddo #else - facvdw=ev1+evdwij - facel=el1+eesij +C MARYSIA + facvdw=(ev1+evdwij)*sss + facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) erij(1)=xj*rmij @@ -3108,9 +3685,9 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 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+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -3149,14 +3726,16 @@ cgrad enddo cgrad enddo do k=1,3 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) + & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) 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) + & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo +C MARYSIA +c endif !sscale 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 @@ -3181,13 +3760,14 @@ C do l=1,2 kkk=kkk+1 muij(kkk)=mu(k,i)*mu(l,j) +c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l #ifdef NEWCORR - gmuij1(kkk)=gtb1(k,i)*mu(l,j) - write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i - gmuij2(kkk)=gUb2(k,i-1)*mu(l,j) - gmuji1(kkk)=mu(k,i)*gtb1(l,j) - write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j - gmuji2(kkk)=mu(k,i)*gUb2(l,j-1) + gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) +c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) + gmuij2(kkk)=gUb2(k,i)*mu(l,j) + gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) +c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) + gmuji2(kkk)=mu(k,i)*gUb2(l,j) #endif enddo enddo @@ -3354,37 +3934,49 @@ cgrad endif C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) +c write (iout,*) 'i',i,' j',j,itype(i),itype(j), +c & ' eel_loc_ij',eel_loc_ij +c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4) C Calculate patrial derivative for theta angle #ifdef NEWCORR geel_loc_ij=a22*gmuij1(1) & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4) - write(iout,*) "derivative over thatai" - write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), - & a33*gmuij1(4) +c write(iout,*) "derivative over thatai" +c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), +c & a33*gmuij1(4) gloc(nphi+i,icg)=gloc(nphi+i,icg)+ & geel_loc_ij*wel_loc - write(iout,*) "derivative over thatai-1" - write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), - & a33*gmuij2(4) - geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3) +c write(iout,*) "derivative over thatai-1" +c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), +c & a33*gmuij2(4) + geel_loc_ij= + & a22*gmuij2(1) + & +a23*gmuij2(2) + & +a32*gmuij2(3) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3) +c Derivative over j residue + geel_loc_ji=a22*gmuji1(1) + & +a23*gmuji1(2) + & +a32*gmuji1(3) & +a33*gmuji1(4) - write(iout,*) "derivative over thataj" - write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), - & a33*gmuji1(4) +c write(iout,*) "derivative over thataj" +c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), +c & a33*gmuji1(4) - gloc(nphi+j,icg)=gloc(nphi+j,icg)+ + gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3) + geel_loc_ji= + & +a22*gmuji2(1) + & +a23*gmuji2(2) + & +a32*gmuji2(3) & +a33*gmuji2(4) - write(iout,*) "derivative over thataj-1" - write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), - & a33*gmuji2(4) +c write(iout,*) "derivative over thataj-1" +c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), +c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc #endif @@ -3392,7 +3984,9 @@ cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij -c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) +c if (eel_loc_ij.ne.0) +c & write (iout,'(a4,2i4,8f9.5)')'chuj', +c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma @@ -3420,14 +4014,14 @@ cgrad enddo cgrad enddo C Remaining derivatives of eello do l=1,3 - gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ - & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4) - gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+ - & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4) - gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+ - & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4) - gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+ - & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4) + gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ + & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) + gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ + & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) + gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ + & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) + gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ + & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) enddo ENDIF C Change 12/26/95 to calculate four-body contributions to H-bonding energy @@ -3640,7 +4234,9 @@ C Third- and fourth-order contributions from turns dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), - & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2) + & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),gpizda1(2,2), + & gpizda2(2,2),auxgmat1(2,2),auxgmatt1(2,2), + & auxgmat2(2,2),auxgmatt2(2,2) double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, @@ -3664,9 +4260,24 @@ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd call checkint_turn3(i,a_temp,eello_turn3_num) call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1)) +c auxalary matices for theta gradient +c auxalary matrix for i+1 and constant i+2 + call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1)) +c auxalary matrix for i+2 and constant i+1 + call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1)) call transpose2(auxmat(1,1),auxmat1(1,1)) + call transpose2(auxgmat1(1,1),auxgmatt1(1,1)) + call transpose2(auxgmat2(1,1),auxgmatt2(1,1)) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) + call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1)) + call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1)) eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) +C Derivatives in theta + gloc(nphi+i,icg)=gloc(nphi+i,icg) + & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 + gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) + & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) cd write (2,*) 'i,',i,' j',j,'eello_turn3', @@ -3740,7 +4351,11 @@ C Third- and fourth-order contributions from turns dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), - & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2) + & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2), + & auxgEvec1(2),auxgEvec2(2),auxgEvec3(2), + & gte1t(2,2),gte2t(2,2),gte3t(2,2), + & gte1a(2,2),gtae3(2,2),gtae3e2(2,2), ae3gte2(2,2), + & gtEpizda1(2,2),gtEpizda2(2,2),gtEpizda3(2,2) double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, @@ -3760,6 +4375,7 @@ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd call checkint_turn4(i,a_temp,eello_turn4_num) c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2 +c write(iout,*)"WCHODZE W PROGRAM" a_temp(1,1)=a22 a_temp(1,2)=a23 a_temp(2,1)=a32 @@ -3771,37 +4387,88 @@ c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) call transpose2(Eug(1,1,i+3),e3t(1,1)) +C Ematrix derivative in theta + call transpose2(gtEUg(1,1,i+1),gte1t(1,1)) + call transpose2(gtEug(1,1,i+2),gte2t(1,1)) + call transpose2(gtEug(1,1,i+3),gte3t(1,1)) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) +c eta1 in derivative theta + call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) +c auxgvec is derivative of Ub2 so i+3 theta call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) +c auxalary matrix of E i+1 + call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1)) c s1=0.0 c gs1=0.0 s1=scalar2(b1(1,i+2),auxvec(1)) -c gs1=scalar2(gtb1(1,i+2),auxgvec(1)) +c derivative of theta i+2 with constant i+3 + gs23=scalar2(gtb1(1,i+2),auxvec(1)) +c derivative of theta i+2 with constant i+2 + gs32=scalar2(b1(1,i+2),auxgvec(1)) +c derivative of E matix in theta of i+1 + gsE13=scalar2(b1(1,i+2),auxgEvec1(1)) + call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) +c ea31 in derivative theta + call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) +c auxilary matrix auxgvec of Ub2 with constant E matirx call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1)) +c auxilary matrix auxgEvec1 of E matix with Ub2 constant + call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1)) + c s2=0.0 c gs2=0.0 s2=scalar2(b1(1,i+1),auxvec(1)) -c gs2=scalar2(gtb1(1,i+1),auxgvec(1)) -c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),ggb1(1,i+2), -c & ggb1(1,i+1) +c derivative of theta i+1 with constant i+3 + gs13=scalar2(gtb1(1,i+1),auxvec(1)) +c derivative of theta i+2 with constant i+1 + gs21=scalar2(b1(1,i+1),auxgvec(1)) +c derivative of theta i+3 with constant i+1 + gsE31=scalar2(b1(1,i+1),auxgEvec3(1)) +c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),gtb1(1,i+2), +c & gtb1(1,i+1) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) +c two derivatives over diffetent matrices +c gtae3e2 is derivative over i+3 + call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1)) +c ae3gte2 is derivative over i+2 + call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) +c three possible derivative over theta E matices +c i+1 + call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1)) +c i+2 + call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1)) +c i+3 + call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) + + gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2)) + gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2)) + gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2)) + eello_turn4=eello_turn4-(s1+s2+s3) +c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)') + & 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3 +cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), +cd & ' eello_turn4_num',8*eello_turn4_num #ifdef NEWCORR -c geel_loc_ij=-(gs1+gs2) -c gloc(nphi+i,icg)=gloc(nphi+i,icg)- -c & gs1 + gloc(nphi+i,icg)=gloc(nphi+i,icg) + & -(gs13+gsE13+gsEE1)*wturn4 + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg) + & -(gs23+gs21+gsEE2)*wturn4 + gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg) + & -(gs32+gsE31+gsEE3)*wturn4 c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- c & gs2 #endif if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn4',i,j,-(s1+s2+s3) -cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), -cd & ' eello_turn4_num',8*eello_turn4_num +c write (iout,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), +c & ' eello_turn4_num',8*eello_turn4_num C Derivatives in gamma(i) call transpose2(EUgder(1,1,i+1),e1tder(1,1)) call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1)) @@ -3969,13 +4636,49 @@ C r0_scp=4.5d0 cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e +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) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) - +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 +c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize +C Condition for being inside the proper box +c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. +c & (xi.lt.((xshift-0.5d0)*boxxsize))) then +c go to 134 +c endif +c 135 continue +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 +c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. +c & (yi.lt.((yshift-0.5d0)*boxysize))) then +c go to 135 +c c endif +c 136 continue +c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize +c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize +cC Condition for being inside the proper box +c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 136 +c endif + 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 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) @@ -3986,10 +4689,75 @@ c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) +c 174 continue +c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize +c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box +c if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 174 +c endif +c 175 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +cC Condition for being inside the proper box +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 175 +c endif +c 176 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +C Condition for being inside the proper box +c if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 176 + 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-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 +C xj=xj-xi +C yj=yj-yi +C zj=zj-zi rij=xj*xj+yj*yj+zj*zj + r0ij=r0_scp r0ijsq=r0ij*r0ij if (rij.lt.r0ijsq) then @@ -4040,6 +4808,9 @@ cgrad enddo enddo ! iint enddo ! i +C enddo !zshift +C enddo !yshift +C enddo !xshift return end C----------------------------------------------------------------------------- @@ -4060,18 +4831,59 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' 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 +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) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) + 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 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 +c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize +C Condition for being inside the proper box +c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. +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 +c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. +c & (yi.lt.((yshift-0.5d0)*boxysize))) then +c go to 135 +c endif +c 136 continue +c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize +c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize +C Condition for being inside the proper box +c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 136 +c endif do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -4082,27 +4894,97 @@ c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + 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 +c 174 continue +c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize +c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box +c if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 174 +c endif +c 175 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +cC Condition for being inside the proper box +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 175 +c endif +c 176 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +C Condition for being inside the proper box +c if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 176 +c endif +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))) fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+e1+e2 + evdw2_14=evdw2_14+(e1+e2)*sss endif evdwij=e1+e2 - evdw2=evdw2+evdwij + evdw2=evdw2+evdwij*sss if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), & bad(itypj,iteli) C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C - fac=-(evdwij+e1)*rrij + fac=-(evdwij+e1)*rrij*sss + fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -4137,10 +5019,14 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - enddo +c endif !endif for sscale cutoff + enddo ! j enddo ! iint enddo ! i +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) @@ -4190,50 +5076,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 @@ -4349,24 +5243,31 @@ c estr=0.0d0 estr1=0.0d0 do i=ibondp_start,ibondp_end - if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then - estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) - do j=1,3 - gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) - & *dc(j,i-1)/vbld(i) - enddo - if (energy_dec) write(iout,*) - & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) - else + if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle +c estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) +c do j=1,3 +c gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) +c & *dc(j,i-1)/vbld(i) +c enddo +c if (energy_dec) write(iout,*) +c & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) +c else +C Checking if it involves dummy (NH3+ or COO-) group + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then +C YES vbldpDUM is the equlibrium length of spring for Dummy atom + diff = vbld(i)-vbldpDUM + else +C NO vbldp0 is the equlibrium lenght of spring for peptide group diff = vbld(i)-vbldp0 - if (energy_dec) write (iout,*) + endif + 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 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) enddo c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) - endif +c endif enddo estr=0.5d0*AKP*estr+estr1 c @@ -4378,7 +5279,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 @@ -4447,7 +5348,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) @@ -4464,7 +5366,7 @@ C Zero the energy function and its derivative at 0 or pi. ichir22=isign(1,itype(i)) endif - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4477,7 +5379,7 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4485,8 +5387,8 @@ C Zero the energy function and its derivative at 0 or pi. z(1)=cos(phii1) #else phii1=phi(i+1) - z(1)=dcos(phii1) #endif + z(1)=dcos(phii1) z(2)=dsin(phii1) else z(1)=0.0D0 @@ -4504,6 +5406,7 @@ C In following comments this theta will be referred to as t_c. bthetk=bthet(k,itype2,ichir21,ichir22) endif thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) +c write(iout,*) 'chuj tu', y(k),z(k) enddo dthett=thet_pred_mean*ssd thet_pred_mean=thet_pred_mean*ss+a0thet(it) @@ -4540,11 +5443,11 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai - if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'ebend',i,ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)') + & 'ebend',i,ethetai,theta(i),itype(i) if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 - gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) + gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo C Ufff.... We've done all this!!! return @@ -4562,7 +5465,8 @@ C--------------------------------------------------------------------------- C Calculate the contributions to both Gaussian lobes. C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time) C The "polynomial part" of the "standard deviation" of this part of -C the distribution. +C the distributioni. +ccc write (iout,*) thetai,thet_pred_mean sig=polthet(3,it) do j=2,0,-1 sig=sig*thet_pred_mean+polthet(j,it) @@ -4592,6 +5496,7 @@ C Following variable is sigma(t_c)**(-2) delthe0=thetai-theta0i term1=-0.5D0*sigcsq*delthec*delthec term2=-0.5D0*sig0inv*delthe0*delthe0 +C write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0 C Following fuzzy logic is to avoid underflows in dexp and subsequent INFs and C NaNs in taking the logarithm. We extract the largest exponent which is added C to the energy (this being the log of the distribution) at the end of energy @@ -4619,6 +5524,7 @@ C Contribution of the bending energy from this theta is just the -log of C the sum of the contributions from the two lobes and the pre-exponential C factor. Simple enough, isn't it? ethetai=(-dlog(termexp)-termm+dlog(termpre)) +C write (iout,*) 'termexp',termexp,termm,termpre,i C NOW the derivatives!!! C 6/6/97 Take into account the deformation. E_theta=(delthec*sigcsq*term1 @@ -4685,7 +5591,11 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle +c print *,i,itype(i-1),itype(i),itype(i-2) + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF + if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -4697,7 +5607,7 @@ C coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4718,7 +5628,7 @@ C propagation of chirality for glycine type sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4850,7 +5760,7 @@ c lprn1=.false. etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg) enddo return end @@ -5616,9 +6526,9 @@ c lprn=.true. do i=iphi_start,iphi_end etors_ii=0.0D0 if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Proline-Proline pair is a special case... @@ -5712,8 +6622,15 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle +C ANY TWO ARE DUMMY ATOMS in row CYCLE +c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. +c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. +c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle + if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF +C For introducing the NH3+ and COO- group please check the etor_d for reference +C and guidance etors_ii=0.0D0 if (iabs(itype(i)).eq.20) then iblock=2 @@ -5814,8 +6731,15 @@ c lprn=.true. etors_d=0.0D0 c write(iout,*) "a tu??" do i=iphid_start,iphid_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle +C ANY TWO ARE DUMMY ATOMS in row CYCLE +C if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. +C & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or. +C & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1)) .or. +C & ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle + if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or. + & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or. + & (itype(i+1).eq.ntyp1)) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5825,9 +6749,21 @@ c write(iout,*) "a tu??" gloci2=0.0D0 iblock=1 if (iabs(itype(i+1)).eq.20) iblock=2 +C Iblock=2 Proline type +C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT +C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO- +C if (itype(i+1).eq.ntyp1) iblock=3 +C The problem of NH3+ group can be resolved by adding new parameters please note if there +C IS or IS NOT need for this +C IF Yes uncomment below and add to parmread.F appropriate changes and to v1cij and so on +C is (itype(i-3).eq.ntyp1) ntblock=2 +C ntblock is N-terminal blocking group C Regular cosine and sine terms do j=1,ntermd_1(itori,itori1,itori2,iblock) +C Example of changes for NH3+ blocking group +C do j=1,ntermd_1(itori,itori1,itori2,iblock,ntblock) +C v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock) v1cij=v1c(1,j,itori,itori1,itori2,iblock) v1sij=v1s(1,j,itori,itori1,itori2,iblock) v2cij=v1c(2,j,itori,itori1,itori2,iblock) @@ -6942,7 +7878,7 @@ C--------------------------------------------------------------------------- if (j.lt.nres-1) then itj1 = itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -7032,14 +7968,14 @@ C parallel orientation of the two CA-CA-CA frames. if (i.gt.1) then iti=itortyp(itype(i)) else - iti=ntortyp+1 + iti=ntortyp endif itk1=itortyp(itype(k+1)) itj=itortyp(itype(j)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif C A1 kernel(j+1) A2T cd do iii=1,2 @@ -7185,7 +8121,7 @@ C Antiparallel orientation of the two CA-CA-CA frames. if (i.gt.1) then iti=itortyp(itype(i)) else - iti=ntortyp+1 + iti=ntortyp endif itk1=itortyp(itype(k+1)) itl=itortyp(itype(l)) @@ -7193,7 +8129,7 @@ C Antiparallel orientation of the two CA-CA-CA frames. if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif C A2 kernel(j-1)T A1T call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), @@ -8151,12 +9087,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, @@ -8327,10 +9263,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 @@ -8347,14 +9283,14 @@ C energy moment and not to the cluster cumulant. if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif itk=itortyp(itype(k)) itk1=itortyp(itype(k+1)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) @@ -8444,7 +9380,7 @@ c---------------------------------------------------------------------------- & auxvec1(2),auxmat1(2,2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C C +C C C Parallel Antiparallel C C C C o o C @@ -8452,10 +9388,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 @@ -8466,19 +9402,19 @@ cd write (2,*) 'eello_graph4: wturn6',wturn6 if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif itk=itortyp(itype(k)) if (k.lt.nres-1) then itk1=itortyp(itype(k+1)) else - itk1=ntortyp+1 + itk1=ntortyp endif itl=itortyp(itype(l)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, @@ -9155,4 +10091,122 @@ 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" +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 + 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=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 + enddo + return + end