X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fenergy_p_new_barrier.F;h=6d6a8179d771190533718b618b6bf9f381dc5001;hb=ae315105cef83bbcab70e2778ef92459690ee784;hp=ddc2a4d30adba53249de9fdcba24166dd3d36ba7;hpb=de4daf245936d67b5c08eb584fca985adab9928e;p=unres.git diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F index ddc2a4d..6d6a817 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -103,10 +103,10 @@ C FG slaves receive the WEIGHTS array wliptran=weights(22) wtube=weights(25) wsaxs=weights(26) - wdfa_dist=weights_(28) - wdfa_tor=weights_(29) - wdfa_nei=weights_(30) - wdfa_beta=weights_(31) + wdfa_dist=weights(28) + wdfa_tor=weights(29) + wdfa_nei=weights(30) + wdfa_beta=weights(31) endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 @@ -118,9 +118,17 @@ c call chainbuild_cart c write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate if (mod(itime_mat,imatupdate).eq.0) then call make_SCp_inter_list +c write (iout,*) "Finished make_SCp_inter_list" +c call flush(iout) call make_SCSC_inter_list +c write (iout,*) "Finished make_SCSC_inter_list" +c call flush(iout) call make_pp_inter_list +c write (iout,*) "Finished make_pp_inter_list" +c call flush(iout) call make_pp_vdw_inter_list +c write (iout,*) "Finished make_pp_vdw_inter_list" +c call flush(iout) endif c print *,'Processor',myrank,' calling etotal ipot=',ipot c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct @@ -169,8 +177,10 @@ C 107 continue #ifdef DFA C BARTEK for dfa test! +c print *,"Processors",MyRank," wdfa",wdfa_dist if (wdfa_dist.gt.0) then call edfad(edfadis) +c print *,"Processors",MyRank," edfadis",edfadis else edfadis=0 endif @@ -366,7 +376,17 @@ c call flush(iout) c write (iout,*) "MULTIBODY_HB ecorr",ecorr,ecorr5,ecorr6,n_corr, c & n_corr1 c call flush(iout) + else + ecorr=0.0d0 + ecorr5=0.0d0 + ecorr6=0.0d0 + eturn6=0.0d0 endif +#else + ecorr=0.0d0 + ecorr5=0.0d0 + ecorr6=0.0d0 + eturn6=0.0d0 #endif c print *,"Processor",myrank," computed Ucorr" c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode @@ -409,15 +429,17 @@ C print *,"za lipidami" call AFMforce(Eafmforce) else if (selfguide.gt.0) then call AFMvel(Eafmforce) + else + Eafmforce=0.0d0 endif if (TUBElog.eq.1) then C print *,"just before call" call calctube(Etube) - elseif (TUBElog.eq.2) then + elseif (TUBElog.eq.2) then call calctube2(Etube) - else - Etube=0.0d0 - endif + else + Etube=0.0d0 + endif #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 @@ -811,7 +833,8 @@ c call flush(iout) #ifdef TIMING c time_allreduce=time_allreduce+MPI_Wtime()-time00 #endif - do i=nnt,nres +c do i=nnt,nres + do i=0,nres do k=1,3 gradbufc(k,i)=0.0d0 enddo @@ -836,7 +859,8 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,-1,-1 +c do i=nres-2,-1,-1 + do i=nres-2,0,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -852,12 +876,13 @@ c enddo #endif #ifdef DEBUG write (iout,*) "gradbufc" - do i=1,nres + do i=0,nres write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) enddo call flush(iout) #endif - do i=-1,nres +c do i=-1,nres + do i=0,nres do j=1,3 gradbufc_sum(j,i)=gradbufc(j,i) gradbufc(j,i)=0.0d0 @@ -866,7 +891,8 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,-1,-1 +c do i=nres-2,-1,-1 + do i=nres-2,0,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -894,7 +920,8 @@ c enddo do k=1,3 gradbufc(k,nres)=0.0d0 enddo - do i=-1,nct +c do i=-1,nct + do i=0,nct do j=1,3 #ifdef SPLITELE C print *,gradbufc(1,13) @@ -999,6 +1026,8 @@ C print *,gradafm(1,13),"AFM" endif #ifdef DEBUG write (iout,*) "gradc gradx gloc after adding" + write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)') + & i,(gradc(j,0,icg),j=1,3),(gradx(j,0,icg),j=1,3) do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)') & i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg) @@ -1028,7 +1057,7 @@ C print *,gradafm(1,13),"AFM" #ifdef MPI if (nfgtasks.gt.1) then do j=1,3 - do i=1,nres + do i=0,nres gradbufc(j,i)=gradc(j,i,icg) gradbufx(j,i)=gradx(j,i,icg) enddo @@ -1055,9 +1084,9 @@ c#undef DEBUG call MPI_Barrier(FG_COMM,IERR) time_barrier_g=time_barrier_g+MPI_Wtime()-time00 time00=MPI_Wtime() - call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres, + call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*(nres+1), & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) - call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres, + call MPI_Reduce(gradbufx(1,0),gradx(1,0,icg),3*(nres+1), & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) @@ -1067,7 +1096,7 @@ c#undef DEBUG time_reduce=time_reduce+MPI_Wtime()-time00 #ifdef DEBUG write (iout,*) "gradc after reduce" - do i=1,nres + do i=0,nres do j=1,3 write (iout,*) i,j,gradc(j,i,icg) enddo @@ -1476,10 +1505,15 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 double precision fcont,fprimcont - double precision sscale,sscagrad + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision gg_lipi(3),gg_lipj(3) double precision boxshift c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -1491,6 +1525,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C Change 12/1/95 num_conti=0 C @@ -1506,6 +1541,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1520,6 +1560,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 + faclip=fac C have you changed here? e1=fac*fac*aa e2=fac*bb @@ -1539,11 +1580,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo cgrad do k=i,j-1 cgrad do l=1,3 @@ -1655,10 +1701,15 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck - double precision sscale,sscagrad double precision boxshift + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision gg_lipi(3),gg_lipj(3) + double precision sscale,sscagrad,sscagradlip,sscalelip c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -1670,6 +1721,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -1681,6 +1733,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1694,6 +1751,7 @@ c do j=istart(i,iint),iend(i,iint) sssgrad1=sscagrad(rij,r_cut_int) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon + faclip=fac C have you changed here? e1=fac*fac*aa e2=fac*bb @@ -1714,11 +1772,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo cgrad do k=i,j-1 cgrad do l=1,3 @@ -1760,13 +1823,16 @@ C integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi, & sss1,sssgrad1 - double precision sscale,sscagrad + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision sscale,sscagrad,sscagradlip,sscalelip double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c if (icall.eq.0) then c lprn=.true. c else @@ -1784,6 +1850,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1822,6 +1889,11 @@ c alf12=0.0D0 yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1844,6 +1916,7 @@ C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives C have you changed here? fac=(rrij*sigsq)**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1871,6 +1944,12 @@ C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1/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 Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad @@ -1911,7 +1990,8 @@ C evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 @@ -1964,7 +2044,8 @@ c write(iout,*) "PO ZWYKLE", evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,' ss' C triple bond artifac removal - do k=j+1,iend(i,iint) +c do k=j+1,iend(i,iint) + do k=j+1,nct C search over all next residues if (dyn_ss_mask(k)) then C check if they are cysteins @@ -2023,9 +2104,15 @@ c alf12=0.0D0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 -C write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) -C if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)') -C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) +c write (iout,*) "aa bb",aa_lip(itypi,itypj), +c & bb_lip(itypi,itypj),aa_aq(itypi,itypj), +c & bb_aq(itypi,itypj),aa,bb +c write (iout,*) (sslipi+sslipj)/2.0d0, +c & (2.0d0-sslipi-sslipj)/2.0d0 + +c write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) +c if (aa.ne.aa_aq(itypi,itypj)) write(iout,'(2e15.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 xj=boxshift(xj-xi,boxxsize) @@ -2053,6 +2140,9 @@ C derivatives. sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij +c if (energy_dec) +c & write (iout,*) "rij",1.0d0/rij," rij_shift",rij_shift, +c & " sig",sig," sig0ij",sig0ij c for diagnostics; uncomment c rij_shift=1.2*sig0ij C I hate to put IF's in the loops, but here don't have another choice!!!! @@ -2061,7 +2151,7 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) - return +c return endif sigder=-sig*sigsq c--------------------------------------------------------------- @@ -2093,8 +2183,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a,2i5,3f10.5)') - & 'r sss evdw',i,j,1.0d0/rij,sss,evdwij + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -2162,7 +2252,8 @@ C double precision dist,sscale,sscagrad,sscagradlip,sscalelip evdw=0.0D0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 @@ -2258,6 +2349,7 @@ 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 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -2285,6 +2377,12 @@ C Calculate gradient components. fac=rij*fac-2*expon*rrij*e_augm fac=fac+(evdwij+e_augm)*sssgrad/sss*rij 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 gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -3450,6 +3548,8 @@ C common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -3556,6 +3656,7 @@ c end if ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -3611,6 +3712,7 @@ c & (zmedi.lt.((-0.5d0)*boxzsize))) then c go to 196 c endif call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -3654,6 +3756,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) C xmedi=xmedi+xshift*boxxsize C ymedi=ymedi+yshift*boxysize C zmedi=zmedi+zshift*boxzsize @@ -3773,6 +3876,9 @@ C------------------------------------------------------------------------------- double precision xmedi,ymedi,zmedi double precision sscale,sscagrad,scalar double precision boxshift + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij, + & faclipij2 + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -3787,6 +3893,7 @@ C 13-go grudnia roku pamietnego... c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j c ind=ind+1 +c write (iout,*) "lipscale",lipscale iteli=itel(i) itelj=itel(j) if (j.eq.i+2 .and. itelj.eq.2) iteli=2 @@ -3807,6 +3914,9 @@ C zj=c(3,j)+0.5D0*dzj-zmedi yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0 + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -3846,14 +3956,15 @@ C fac_shield(j)=0.6 el1=el1*fac_shield(i)**2*fac_shield(j)**2 el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*sss*faclipij2 else fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) - ees=ees+eesij*sss + ees=ees+eesij*sss*faclipij2 endif - evdw1=evdw1+evdwij*sss + ees=ees + evdw1=evdw1+evdwij*sss*faclipij2 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, @@ -3862,8 +3973,9 @@ cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,3f10.5)') & 'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss,rij - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - & fac_shield(i),fac_shield(j) + write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij, + & faclipij2 endif C @@ -3880,7 +3992,7 @@ C * * Radial derivatives. First process both termini of the fragment (i,j) * - aux=facel*sss+rmij*sssgrad*eesij + aux=(facel*sss+rmij*sssgrad*eesij)*faclipij2 ggg(1)=aux*xj ggg(2)=aux*yj ggg(3)=aux*zj @@ -3962,15 +4074,14 @@ c 9/28/08 AL Gradient compotents will be summed only at the end C print *,"before", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc_long(k,j)=gelc_long(k,j)+ggg(k) -C & +grad_shield(k,j)*eesij/fac_shield(j) gelc_long(k,i)=gelc_long(k,i)-ggg(k) -C & +grad_shield(k,i)*eesij/fac_shield(i) -C gelc_long(k,i-1)=gelc_long(k,i-1) -C & +grad_shield(k,i)*eesij/fac_shield(i) -C gelc_long(k,j-1)=gelc_long(k,j-1) -C & +grad_shield(k,j)*eesij/fac_shield(j) enddo -C print *,"bafter", gelc_long(1,i), gelc_long(1,j) + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2*sss + + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2*sss + * * Loop over residues i+1 thru j-1. @@ -3980,7 +4091,7 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - facvdw=facvdw+sssgrad*rmij*evdwij + facvdw=(facvdw+sssgrad*rmij*evdwij)*faclipij2 ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj @@ -3994,6 +4105,11 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 * * Loop over residues i+1 thru j-1. * @@ -4004,7 +4120,7 @@ cgrad enddo cgrad enddo #else C MARYSIA - facvdw=(ev1+evdwij) + facvdw=(ev1+evdwij)*faclipij2 facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel)*sss @@ -4047,6 +4163,10 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 #endif * * Angular part @@ -4064,7 +4184,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* - & fac_shield(i)**2*fac_shield(j)**2*sss + & fac_shield(i)**2*fac_shield(j)**2*sss*faclipij2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -4085,11 +4205,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -4324,7 +4444,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') c & 'eelloc',i,j,eel_loc_ij C Now derivative over eel_loc @@ -4382,7 +4502,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -4398,7 +4518,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -4411,7 +4531,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -4423,7 +4543,7 @@ 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 - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -4439,12 +4559,12 @@ C Partial derivatives in virtual-bond dihedral angles gamma & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) aux=eel_loc_ij/sss*sssgrad*rmij ggg(1)=aux*xj @@ -4453,13 +4573,19 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij + + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij + cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -4469,19 +4595,19 @@ 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)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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)) - & *fac_shield(i)*fac_shield(j)*sss + & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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)) - & *fac_shield(i)*fac_shield(j)*sss + & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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)) - & *fac_shield(i)*fac_shield(j)*sss + & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) + & *fac_shield(i)*fac_shield(j)*sss*faclipij enddo ENDIF @@ -4734,6 +4860,8 @@ C Third- and fourth-order contributions from turns common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij j=i+2 c write (iout,*) "eturn3",i,j,j1,j2 a_temp(1,1)=a22 @@ -4771,7 +4899,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) if (energy_dec) write (iout,'(6heturn3,2i5,0pf7.3)') i,i+2, @@ -4780,10 +4908,10 @@ C#ifdef NEWCORR C Derivatives in theta gloc(nphi+i,icg)=gloc(nphi+i,icg) & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C#endif C Derivatives in shield mode @@ -4838,14 +4966,14 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Cartesian derivatives do l=1,3 c ghalf1=0.5d0*agg(l,1) @@ -4859,7 +4987,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) @@ -4868,7 +4996,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -4876,7 +5004,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -4884,8 +5012,17 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo + gshieldc_t3(3,i)=gshieldc_t3(3,i)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j)=gshieldc_t3(3,j)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + return end C------------------------------------------------------------------------------- @@ -5014,7 +5151,7 @@ C fac_shield(i)=0.6 C fac_shield(j)=0.4 endif eello_turn4=eello_turn4-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) @@ -5062,12 +5199,6 @@ C & *2.0 & grad_shield(k,j)*eello_t4/fac_shield(j) enddo endif - - - - - - cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), cd & ' eello_turn4_num',8*eello_turn4_num #ifdef NEWCORR @@ -5097,7 +5228,7 @@ C Derivatives in gamma(i) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -5106,7 +5237,7 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -5118,7 +5249,7 @@ C Derivatives in gamma(i+2) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -5138,7 +5269,7 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif C Remaining derivatives of this turn contribution @@ -5157,7 +5288,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -5172,7 +5303,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -5187,7 +5318,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -5203,8 +5334,16 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) - enddo + & *fac_shield(i)*fac_shield(j)*faclipij + enddo + gshieldc_t4(3,i)=gshieldc_t4(3,i)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j)=gshieldc_t4(3,j)+ + & ssgradlipj*eello_t4/4.0d0*lipscale + gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+ + & ssgradlipj*eello_t4/4.0d0*lipscale return end C----------------------------------------------------------------------------- @@ -7691,16 +7830,20 @@ c enddo c min_odl=minval(distancek) - do kk=1,constr_homology - if(l_homo(kk,ii)) then - min_odl=distancek(kk) - exit - endif - enddo - do kk=1,constr_homology - if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + if (nexl.gt.0) then + min_odl=0.0d0 + else + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) & min_odl=distancek(kk) - enddo + enddo + endif c write (iout,* )"min_odl",min_odl #ifdef DEBUG @@ -11778,6 +11921,7 @@ C--bufliptop--- here true lipid starts C lipid C--buflipbot--- lipid ends buffore starts C--bordlipbot--buffore ends +c call cartprint eliptran=0.0 do i=ilip_start,ilip_end C do i=1,1 @@ -11832,6 +11976,8 @@ CV do i=1,1 if (itype(i).eq.ntyp1) cycle positi=(mod(c(3,i+nres),boxzsize)) if (positi.le.0) positi=positi+boxzsize +c write(iout,*) "i",i," positi",positi,bordlipbot,buflipbot, +c & bordliptop 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) @@ -11842,6 +11988,8 @@ C the energy transfer exist if (positi.lt.buflipbot) then fracinbuf=1.0d0- & ((positi-bordlipbot)/lipbufthick) +c write (iout,*) "i",i,itype(i)," fracinbuf",fracinbuf +c write (iout,*) "i",i," liptranene",liptranene(itype(i)) C lipbufthick is thickenes of lipid buffore sslip=sscalelip(fracinbuf) ssgradlip=-sscagradlip(fracinbuf)/lipbufthick @@ -13217,11 +13365,16 @@ c-------------------------------------------------------------------------- subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi) implicit none include 'DIMENSIONS' + include 'COMMON.IOUNITS' include 'COMMON.CHAIN' double precision xi,yi,zi,sslipi,ssgradlipi double precision fracinbuf double precision sscalelip,sscagradlip - +#ifdef DEBUG + write (iout,*) "bordlipbot",bordlipbot," bordliptop",bordliptop + write (iout,*) "buflipbot",buflipbot," lipbufthick",lipbufthick + write (iout,*) "xi yi zi",xi,yi,zi +#endif if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then C the energy transfer exist if (zi.lt.buflipbot) then @@ -13242,5 +13395,8 @@ C lipbufthick is thickenes of lipid buffore sslipi=0.0d0 ssgradlipi=0.0 endif +#ifdef DEBUG + write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi +#endif return end