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=9ec8107a38505730ffcacf204a23584d50c4c37e;hb=7a9066ef8247b567cb33f4f6f3483e076614c7c4;hp=b36b9a85b829039e0f8cead26aeefdbe3a1a96d8;hpb=c711143ad3fffb04d27b55aa823f399b8343c4c5;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 b36b9a8..9ec8107 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -30,12 +30,14 @@ c include 'COMMON.MD' include 'COMMON.SPLITELE' include 'COMMON.TORCNSTR' include 'COMMON.SAXS' + include 'COMMON.MD' double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc, & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr, & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6, & eliptran,Eafmforce,Etube, & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet integer n_corr,n_corr1 + double precision time01 #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -66,7 +68,8 @@ C FG slaves as WEIGHTS array. weights_(17)=wbond weights_(18)=scal14 weights_(21)=wsccor - weights_(22)=wtube + weights_(22)=wliptran + weights_(25)=wtube weights_(26)=wsaxs weights_(28)=wdfa_dist weights_(29)=wdfa_tor @@ -98,23 +101,42 @@ C FG slaves receive the WEIGHTS array wbond=weights(17) scal14=weights(18) wsccor=weights(21) - wtube=weights(22) + 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 c call chainbuild_cart endif -#ifndef DFA - edfadis=0.0d0 - edfator=0.0d0 - edfanei=0.0d0 - edfabet=0.0d0 + if (nfgtasks.gt.1) then + call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR) + endif +c write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate + if (mod(itime_mat,imatupdate).eq.0) then +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif + 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) +c call make_pp_vdw_inter_list +c write (iout,*) "Finished make_pp_vdw_inter_list" +c call flush(iout) +#ifdef TIMING_ENE + time_list=time_list+MPI_Wtime()-time01 #endif + endif c print *,'Processor',myrank,' calling etotal ipot=',ipot c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct #else @@ -125,10 +147,20 @@ c endif #ifdef TIMING time00=MPI_Wtime() #endif + +#ifndef DFA + edfadis=0.0d0 + edfator=0.0d0 + edfanei=0.0d0 + edfabet=0.0d0 +#endif C C Compute the side-chain and electrostatic interaction energy C C print *,ipot +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj(evdw) @@ -153,10 +185,15 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +#ifdef TIMING_ENE + time_evdw=time_evdw+MPI_Wtime()-time01 +#endif #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 @@ -192,6 +229,9 @@ c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time_vec=time_vec+MPI_Wtime()-time01 #endif +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif C Introduction of shielding effect first for each peptide group C the shielding factor is set this factor is describing how each C peptide group is shielded by side-chains @@ -228,6 +268,9 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, c & eello_turn4) endif +#ifdef TIMING_ENE + time_eelec=time_eelec+MPI_Wtime()-time01 +#endif c#ifdef TIMING c time_enecalc=time_enecalc+MPI_Wtime()-time00 c#endif @@ -236,6 +279,9 @@ C C Calculate excluded-volume interaction energy between peptide groups C and side chains. C +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif if (ipot.lt.6) then if(wscp.gt.0d0) then call escp(evdw2,evdw2_14) @@ -247,6 +293,9 @@ C c write (iout,*) "Soft-sphere SCP potential" call escp_soft_sphere(evdw2,evdw2_14) endif +#ifdef TIMING_ENE + time_escp=time_escp+MPI_Wtime()-time01 +#endif c c Calculate the bond-stretching energy c @@ -352,7 +401,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 @@ -387,21 +446,25 @@ C based on partition function C print *,"przed lipidami" if (wliptran.gt.0) then call Eliptransfer(eliptran) + else + eliptran=0.0d0 endif C print *,"za lipidami" if (AFMlog.gt.0) then 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 @@ -795,17 +858,18 @@ 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 enddo -#ifdef DEBUG - write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end - write (iout,*) (i," jgrad_start",jgrad_start(i), - & " jgrad_end ",jgrad_end(i), - & i=igrad_start,igrad_end) -#endif +c#ifdef DEBUG +c write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end +c write (iout,*) (i," jgrad_start",jgrad_start(i), +c & " jgrad_end ",jgrad_end(i), +c & i=igrad_start,igrad_end) +c#endif c c Obsolete and inefficient code; we can make the effort O(n) and, therefore, c do not parallelize this part. @@ -820,7 +884,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 @@ -836,12 +901,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 @@ -850,7 +916,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 @@ -878,7 +945,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) @@ -983,6 +1051,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) @@ -1012,7 +1082,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 @@ -1039,9 +1109,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) @@ -1051,7 +1121,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 @@ -1449,45 +1519,74 @@ C include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' + include 'COMMON.SPLITELE' #ifdef FOURBODY include 'COMMON.CONTACTS' include 'COMMON.CONTMAT' #endif double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,num_conti,iint + integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, - & sigij,r0ij,rcut + & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 double precision fcont,fprimcont + 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 + external boxshift c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 - do i=iatsc_s,iatsc_e + 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) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) 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 C Calculate SC interaction energy. C - do iint=1,nint_gr(i) +c do iint=1,nint_gr(i) cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) - do j=istart(i,iint),iend(i,iint) +c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle - 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) + 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) C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj rrij=1.0D0/rij + sqrij=dsqrt(rij) + sss1=sscale(sqrij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(sqrij,r_cut_int) + 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 @@ -1498,19 +1597,25 @@ cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') 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 + evdw=evdw+sss1*evdwij C C Calculate the components of the gradient in DC and X C - fac=-rrij*(e1+evdwij) + fac=-rrij*(e1+evdwij)*sss1 + & +evdwij*sssgrad1/sqrij/expon 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 @@ -1575,8 +1680,8 @@ cd & i,j,(gacont(kk,num_conti,i),kk=1,3) endif endif #endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint C Change 12/1/95 #ifdef FOURBODY num_cont(i)=num_conti @@ -1615,38 +1720,64 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' + include 'COMMON.SPLITELE' double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,iint + integer i,j,k,itypi,itypj,itypi1,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, - & fac_augm,e_augm,r_inv_ij,r_shift_inv + & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck + 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 - do i=iatsc_s,iatsc_e + 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) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) 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 - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle - 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) + 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) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij + sss1=sscale(rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + 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 @@ -1658,27 +1789,33 @@ cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - evdw=evdw+evdwij + evdw=evdw+evdwij*sss1 C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) + & +evdwij*sssgrad1*r_inv_ij/expon 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 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) cgrad enddo cgrad enddo - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i do i=1,nct do j=1,3 @@ -1705,29 +1842,41 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SPLITELE' integer icall common /srutu/ icall double precision evdw - integer itypi,itypj,itypi1,iint,ind - double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi + integer itypi,itypj,itypi1,iint,ind,ikont + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi, + & sss1,sssgrad1 + 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 lprn=.false. c endif ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) 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) @@ -1736,8 +1885,8 @@ c dsci_inv=dsc_inv(itypi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1762,9 +1911,18 @@ 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) + 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) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1775,19 +1933,23 @@ cd else cd rrij=rrsave(ind) cd endif rij=dsqrt(rrij) + sss1=sscale(1.0d0/rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(1.0d0/rij,r_cut_int) 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 + faclip=fac 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 + evdw=evdw+sss1*evdwij if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa @@ -1803,15 +1965,22 @@ C Calculate gradient components. fac=-expon*(e1+evdwij) sigder=fac/sigsq fac=rrij*fac + & +evdwij*sssgrad1/sss1*rij 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 - enddo ! j - enddo ! iint +! enddo ! j +! enddo ! iint enddo ! i c stop return @@ -1837,18 +2006,23 @@ C include 'COMMON.SPLITELE' include 'COMMON.SBRIDGE' logical lprn - integer xshift,yshift,zshift,subchap double precision evdw - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, - & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, - & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip + double precision boxshift + double precision x0,y0,r012,rij12,facx0, + & facx02,afacx0,bfacx0,abfacx0,Afac,BBfac,Afacsig,Bfacsig +c alpha_GB=0.5d0 +c alpha_GB=0.01d0 +c alpha_GB1=1.0d0+1.0d0/alpha_GB 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 @@ -1857,73 +2031,24 @@ 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 +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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 + call to_box(xi,yi,zi) 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 - + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C xi=xi+xshift*boxxsize C yi=yi+yshift*boxysize C zi=zi+zshift*boxzsize @@ -1938,27 +2063,29 @@ c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN c write(iout,*) "PRZED ZWYKLE", evdwij call dyn_ssbond_ene(i,j,evdwij) c write(iout,*) "PO ZWYKLE", evdwij +c call flush(iout) evdw=evdw+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 + if (dyn_ss_mask(k)) then C check if they are cysteins C write(iout,*) 'k=',k c write(iout,*) "PRZED TRI", evdwij - evdwij_przed_tri=evdwij - call triple_ssbond_ene(i,j,k,evdwij) + evdwij_przed_tri=evdwij + call triple_ssbond_ene(i,j,k,evdwij) c if(evdwij_przed_tri.ne.evdwij) then c write (iout,*) "TRI:", evdwij, evdwij_przed_tri c endif @@ -1966,30 +2093,30 @@ c endif c write(iout,*) "PO TRI", evdwij C call the energy function that removes the artifical triple disulfide C bond the soubroutine is located in ssMD.F - evdw=evdw+evdwij - if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,'tss' - endif!dyn_ss_mask(k) - enddo! k + endif!dyn_ss_mask(k) + enddo! k ELSE - ind=ind+1 - itypj=iabs(itype(j)) - if (itypj.eq.ntyp1) cycle + ind=ind+1 + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) - dscj_inv=vbld_inv(j+nres) + dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, c & 1.0d0/vbld(j+nres) c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) - sig0ij=sigma(itypi,itypj) - chi1=chi(itypi,itypj) - chi2=chi(itypj,itypi) - chi12=chi1*chi2 - chip1=chip(itypi) - chip2=chip(itypj) - chip12=chip1*chip2 - alf1=alp(itypi) - alf2=alp(itypj) - alf12=0.5D0*(alf1+alf2) + sig0ij=sigma(itypi,itypj) + chi1=chi(itypi,itypj) + chi2=chi(itypj,itypi) + chi12=chi1*chi2 + chip1=chip(itypi) + chip2=chip(itypj) + chip12=chip1*chip2 + alf1=alp(itypi) + alf2=alp(itypj) + alf12=0.5D0*(alf1+alf2) C For diagnostics only!!! c chi1=0.0D0 c chi2=0.0D0 @@ -2000,195 +2127,207 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - 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 (zj.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) - sslipj=sscalelip(fracinbuf) - ssgradlipj=sscagradlip(fracinbuf)/lipbufthick - else - sslipj=1.0d0 - ssgradlipj=0.0 - endif - else - sslipj=0.0d0 - ssgradlipj=0.0 - endif - aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 - bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 -C write(iout,*) "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)) + xj=c(1,nres+j) + 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 +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 - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 - xj=xj_safe+xshift*boxxsize - yj=yj_safe+yshift*boxysize - zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - subchap=1 - endif - enddo - enddo - enddo - if (subchap.eq.1) then - xj=xj_temp-xi - yj=yj_temp-yi - zj=zj_temp-zi - else - xj=xj_safe-xi - yj=yj_safe-yi - zj=zj_safe-zi - endif - dxj=dc_norm(1,nres+j) - dyj=dc_norm(2,nres+j) - dzj=dc_norm(3,nres+j) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) + 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)) - + rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + rij=dsqrt(rrij) + sss=sscale(1.0d0/rij,r_cut_int) c write (iout,'(a7,4f8.3)') c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb - if (sss.gt.0.0d0) then + if (sss.eq.0.0d0) cycle + sssgrad=sscagrad(1.0d0/rij,r_cut_int) C Calculate angle-dependent terms of energy and contributions to their C derivatives. - call sc_angular - sigsq=1.0D0/sigsq - sig=sig0ij*dsqrt(sigsq) - rij_shift=1.0D0/rij-sig+sig0ij + call sc_angular + 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!!!! - if (rij_shift.le.0.0D0) then - evdw=1.0D20 +c if (rij_shift.le.0.0D0) then + x0=alpha_GB*(sig-sig0ij) + if (energy_dec) write (iout,*) i,j," x0",x0 + if (rij_shift.le.x0) then +c sig=2.0d0*sig0ij + sigder=-sig*sigsq +c sigder=0.0d0 + fac=rij**expon + rij12=fac*fac +c rij12=1.0d0 + x0=alpha_GB*(sig-sig0ij) + facx0=1.0d0/x0**expon + facx02=facx0*facx0 + r012=((1.0d0+alpha_GB)*(sig-sig0ij))**(2*expon) + afacx0=aa*facx02 + bfacx0=bb*facx0 + abfacx0=afacx0+0.5d0*bfacx0 + Afac=alpha_GB1*abfacx0 + Afacsig=0.5d0*alpha_GB1*bfacx0/(sig-sig0ij) + BBfac=Afac-(afacx0+bfacx0) +c BBfac=0.0d0 + Bfacsig=(-alpha_GB1*(abfacx0+afacx0)+ + & (afacx0+afacx0+bfacx0))/(sig-sig0ij) +c Bfacsig=0.0d0 + Afac=Afac*r012 + Afacsig=Afacsig*r012 +c Afac=1.0d0 +c Afacsig=0.0d0 +c w(x)=4*eps*((1.0+1.0/alpha_GB)*(y0**12-0.5*y0**6)*(r0/x)**12-(1+1/alpha)*(y0**12-0.5*y0**6)+y0**12-y0**6) +c eps1=1.0d0 +c eps2rt=1.0d0 +c eps3rt=1.0d0 + e1 = eps1*eps2rt*eps3rt*Afac*rij12 + e2 = -eps1*eps2rt*eps3rt*BBfac + evdwij = e1+e2 + eps2der=evdwij*eps3rt + eps3der=evdwij*eps2rt +c eps2der=0.0d0 +c eps3der=0.0d0 +c eps1_om12=0.0d0 + evdwij=evdwij*eps2rt*eps3rt +c Afacsig=0.d0 +c Bfacsig=0.0d0 + if (lprn) then + write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa + write (iout,'(2(a3,i3,2x),18(0pf9.5))') + & restyp(itypi),i,restyp(itypj),j, + & epsi,sigm,chi1,chi2,chip1,chip2, + & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, + & eps1*eps2rt**2*eps3rt**2,om1,om2,om12, + & 1.0D0/rij,rij_shift, + & evdwij + endif + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'RE r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij + evdw=evdw+evdwij*sss +C Calculate gradient components. + e1=e1*eps2rt*eps3rt + sigder=-expon*eps1*eps2rt*eps2rt*eps3rt*eps3rt + & *(Afacsig*rij12-Bfacsig)*sigder + fac=-2.0d0*expon*e1*rij*rij +c print '(2i4,6f8.4)',i,j,sss,sssgrad* +c & evdwij,fac,sigma(itypi,itypj),expon + fac=fac+evdwij*sssgrad/sss*rij +c fac=0.0d0 +c write (iout,*) "sigder",sigder," fac",fac," e1",e1, +c & " e2",e2," sss",sss," sssgrad",sssgrad,"esp123", +c & eps1*eps2rt**2*eps3rt**2 +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 evdw=1.0D20 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 - endif - sigder=-sig*sigsq +c return + else + rij_shift=1.0D0/rij_shift + sigder=-sig*sigsq c--------------------------------------------------------------- - rij_shift=1.0D0/rij_shift - fac=rij_shift**expon + fac=rij_shift**expon 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 + faclip=fac + e1=fac*fac*aa + e2=fac*bb + evdwij=eps1*eps2rt*eps3rt*(e1+e2) + eps2der=evdwij*eps3rt + eps3der=evdwij*eps2rt C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij, C &((sslipi+sslipj)/2.0d0+ C &(2.0d0-sslipi-sslipj)/2.0d0) c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 - evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij*sss - if (lprn) then - sigm=dabs(aa/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, - & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, - & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, - & evdwij - endif + evdwij=evdwij*eps2rt*eps3rt + evdw=evdw+evdwij*sss + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'GB r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij + if (lprn) then + write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij + 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, + & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, + & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, + & evdwij + endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij C Calculate gradient components. - e1=e1*eps1*eps2rt**2*eps3rt**2 - 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 + e1=e1*eps1*eps2rt**2*eps3rt**2 + 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*sssgrad/sss*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 + 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 + endif C Calculate angular part of the gradient. - call sc_grad - endif +c call sc_grad_scale(sss) + call sc_grad ENDIF ! dyn_ss - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i C enddo ! zshift C enddo ! yshift @@ -2214,66 +2353,43 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' - integer xshift,yshift,zshift,subchap + include 'COMMON.SPLITELE' + double precision boxshift integer icall common /srutu/ icall logical lprn double precision evdw - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, & xi,yi,zi,fac_augm,e_augm double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, - & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, - & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip,sssgrad1 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 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) 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 + call to_box(xi,yi,zi) 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 - + 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) @@ -2282,8 +2398,8 @@ c dsci_inv=dsc_inv(itypi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -2310,129 +2426,86 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 -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- - & ((zj-bordlipbot)/lipbufthick) -C lipbufthick is thickenes of lipid buffore - sslipj=sscalelip(fracinbuf) - ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zj.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) - sslipj=sscalelip(fracinbuf) - ssgradlipj=sscagradlip(fracinbuf)/lipbufthick - else - sslipj=1.0d0 - ssgradlipj=0.0 - endif - else - sslipj=0.0d0 - ssgradlipj=0.0 - endif - aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 - bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 + xj=c(1,nres+j) + 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 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,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj - 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) - rrij=1.0D0/(xj*xj+yj*yj+zj*zj) - rij=dsqrt(rrij) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) + dxj=dc_norm(1,nres+j) + dyj=dc_norm(2,nres+j) + dzj=dc_norm(3,nres+j) + rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + rij=dsqrt(rrij) + sss=sscale(1.0d0/rij,r_cut_int) + if (sss.eq.0.0d0) cycle + sssgrad=sscagrad(1.0d0/rij,r_cut_int) C Calculate angle-dependent terms of energy and contributions to their C derivatives. - call sc_angular - sigsq=1.0D0/sigsq - sig=sig0ij*dsqrt(sigsq) - rij_shift=1.0D0/rij-sig+r0ij + call sc_angular + sigsq=1.0D0/sigsq + sig=sig0ij*dsqrt(sigsq) + rij_shift=1.0D0/rij-sig+r0ij C I hate to put IF's in the loops, but here don't have another choice!!!! - if (rij_shift.le.0.0D0) then - evdw=1.0D20 - return - endif - sigder=-sig*sigsq + if (rij_shift.le.0.0D0) then + evdw=1.0D20 + return + endif + sigder=-sig*sigsq c--------------------------------------------------------------- - rij_shift=1.0D0/rij_shift - fac=rij_shift**expon - e1=fac*fac*aa - e2=fac*bb - evdwij=eps1*eps2rt*eps3rt*(e1+e2) - eps2der=evdwij*eps3rt - eps3der=evdwij*eps2rt - fac_augm=rrij**expon - e_augm=augm(itypi,itypj)*fac_augm - evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij+e_augm - if (lprn) then - sigm=dabs(aa/bb)**(1.0D0/6.0D0) - epsi=bb**2/aa - write (iout,'(2(a3,i3,2x),17(0pf7.3))') + 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) + eps2der=evdwij*eps3rt + eps3der=evdwij*eps2rt + fac_augm=rrij**expon + e_augm=augm(itypi,itypj)*fac_augm + evdwij=evdwij*eps2rt*eps3rt + evdw=evdw+evdwij+e_augm + if (lprn) then + 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), & chi1,chi2,chip1,chip2, & eps1,eps2rt**2,eps3rt**2, & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, & evdwij+e_augm - endif + endif C Calculate gradient components. - e1=e1*eps1*eps2rt**2*eps3rt**2 - 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 + e1=e1*eps1*eps2rt**2*eps3rt**2 + fac=-expon*(e1+evdwij)*rij_shift + sigder=fac*sigder + 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(1)=xj*fac - gg(2)=yj*fac - gg(3)=zj*fac + 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 C Calculate angular part of the gradient. - call sc_grad - enddo ! j - enddo ! iint +c call sc_grad_scale(sss) + call sc_grad +c enddo ! j +c enddo ! iint enddo ! i end C----------------------------------------------------------------------------- @@ -2581,27 +2654,32 @@ C include 'COMMON.IOUNITS' c include 'COMMON.CONTACTS' dimension gg(3) + double precision boxshift cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct evdw=0.0D0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do ikont=g_listscsc_start,g_listscsc_end + i=newcontlisti(ikont) + j=newcontlistj(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + call to_box(xi,yi,zi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) +c do iint=1,nint_gr(i) cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) - do j=istart(i,iint),iend(i,iint) +c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi + xj=boxshift(c(1,nres+j)-xi,boxxsize) + yj=boxshift(c(2,nres+j)-yi,boxysize) + zj=boxshift(c(3,nres+j)-zi,boxzsize) rij=xj*xj+yj*yj+zj*zj c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj r0ij=r0(itypi,itypj) @@ -2632,8 +2710,8 @@ cgrad do l=1,3 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) cgrad enddo cgrad enddo - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i return end @@ -2658,7 +2736,7 @@ c include 'COMMON.CONTACTS' include 'COMMON.VECTORS' include 'COMMON.FFIELD' dimension ggg(3) - integer xshift,yshift,zshift + double precision boxshift C write(iout,*) 'In EELEC_soft_sphere' ees=0.0D0 evdw1=0.0D0 @@ -2674,12 +2752,7 @@ C 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 + call to_box(xmedi,ymedi,zmedi) num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -2696,46 +2769,13 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) 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 + call to_box(xj,yj,zj) + xj=boxshift(xj-xmedi,boxxsize) + yj=boxshift(yj-ymedi,boxysize) + zj=boxshift(zj-zmedi,boxzsize) rij=xj*xj+yj*yj+zj*zj - sss=sscale(sqrt(rij)) - sssgrad=sscagrad(sqrt(rij)) + sss=sscale(sqrt(rij),r_cut_int) + sssgrad=sscagrad(sqrt(rij),r_cut_int) if (rij.lt.r0ijsq) then evdw1ij=0.25d0*(rij-r0ijsq)**2 fac=rij-r0ijsq @@ -3124,7 +3164,7 @@ c write(iout,*) 'b2=',(b2(k,i-2),k=1,2) #endif enddo - mu=0.0d0 + mu(:,:nres)=0.0d0 #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else @@ -3623,6 +3663,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/ @@ -3728,12 +3770,8 @@ c end if 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 + 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) @@ -3788,13 +3826,8 @@ 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 - + call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -3814,7 +3847,10 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c CTU KURWA - do i=iatel_s,iatel_e +c do i=iatel_s,iatel_e + do ikont=g_listpp_start,g_listpp_end + i=newcontlistppi(ikont) + j=newcontlistppj(ikont) C do i=75,75 c if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 @@ -3834,12 +3870,8 @@ c & .or. itype(i-1).eq.ntyp1 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 + 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 @@ -3875,11 +3907,11 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) #endif C I TU KURWA - do j=ielstart(i),ielend(i) +c do j=ielstart(i),ielend(i) C do j=16,17 C write (iout,*) i,j C if (j.le.1) cycle - if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 + if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds c & .or.((j+2).gt.nres) c & .or.((j-1).le.0) @@ -3887,8 +3919,8 @@ C end of changes by Ana c & .or.itype(j+2).eq.ntyp1 c & .or.itype(j-1).eq.ntyp1 &) cycle - call eelecij(i,j,ees,evdw1,eel_loc) - enddo ! j + call eelecij(i,j,ees,evdw1,eel_loc) +c enddo ! j #ifdef FOURBODY num_cont_hb(i)=num_conti #endif @@ -3909,7 +3941,7 @@ cd print *,"Processor",fg_rank," t_eelecij",t_eelecij end C------------------------------------------------------------------------------- subroutine eelecij(i,j,ees,evdw1,eel_loc) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" @@ -3933,14 +3965,35 @@ C------------------------------------------------------------------------------- include 'COMMON.TIME1' include 'COMMON.SPLITELE' include 'COMMON.SHIELD' - dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), + double precision 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), & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4), & gmuij2(4),gmuji2(4) + double precision dxi,dyi,dzi + double precision dx_normi,dy_normi,dz_normi,aux + integer j1,j2,lll,num_conti 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 + integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ilist,iresshield + double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp + double precision ees,evdw1,eel_loc,aaa,bbb,ael3i + double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj, + & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4, + & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa, + & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der, + & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij, + & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp, + & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp, + & ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield + double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji + 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/ @@ -3952,10 +4005,10 @@ C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ - integer xshift,yshift,zshift 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 @@ -3975,77 +4028,20 @@ 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 + 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) 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)) + sss=sscale(dsqrt(rij),r_cut_int) + if (sss.eq.0.0d0) return + sssgrad=sscagrad(dsqrt(rij),r_cut_int) c if (sss.gt.0.0d0) then rrmij=1.0D0/rij rij=dsqrt(rij) @@ -4075,25 +4071,26 @@ 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 + 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, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') - &'evdw1',i,j,evdwij - &,iteli,itelj,aaa,evdw1,sss - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - &fac_shield(i),fac_shield(j) + 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,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij, + & faclipij2 endif C @@ -4110,9 +4107,10 @@ C * * Radial derivatives. First process both termini of the fragment (i,j) * - ggg(1)=facel*xj - ggg(2)=facel*yj - ggg(3)=facel*zj + aux=(facel*sss+rmij*sssgrad*eesij)*faclipij2 + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j @@ -4146,10 +4144,10 @@ C endif iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) - & *2.0 + & *2.0*sss gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) @@ -4172,13 +4170,13 @@ C endif do k=1,3 gshieldc(k,i)=gshieldc(k,i)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss gshieldc(k,j)=gshieldc(k,j)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss gshieldc(k,i-1)=gshieldc(k,i-1)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss gshieldc(k,j-1)=gshieldc(k,j-1)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss enddo endif @@ -4191,15 +4189,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. @@ -4209,15 +4206,10 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - 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 + facvdw=(facvdw+sssgrad*rmij*evdwij)*faclipij2 + ggg(1)=facvdw*xj + ggg(2)=facvdw*yj + ggg(3)=facvdw*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -4228,6 +4220,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. * @@ -4238,10 +4235,11 @@ cgrad enddo cgrad enddo #else C MARYSIA - facvdw=(ev1+evdwij)*sss + facvdw=(ev1+evdwij)*faclipij2 facel=(el1+eesij) fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel) + fac=-3*rrmij*(facvdw+facvdw+facel)*sss + & +(evdwij+eesij)*sssgrad*rrmij erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij @@ -4280,6 +4278,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 @@ -4297,7 +4299,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 + & fac_shield(i)**2*fac_shield(j)**2*sss*faclipij2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -4317,12 +4319,12 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) 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)) - & *fac_shield(i)**2*fac_shield(j)**2 + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss + & *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)) - & *fac_shield(i)**2*fac_shield(j)**2 + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss + & *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 @@ -4557,7 +4559,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) + & *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 @@ -4615,7 +4617,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j) + & *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) @@ -4631,7 +4633,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) + & *fac_shield(i)*fac_shield(j)*sss*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -4644,7 +4646,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -4656,7 +4658,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) + & *fac_shield(i)*fac_shield(j)*sss*faclipij #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -4672,23 +4674,33 @@ 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) + & *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) + & *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 + ggg(2)=aux*yj + ggg(3)=aux*zj do l=1,3 - ggg(l)=(agg(l,1)*muij(1)+ + 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) + & *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) @@ -4698,19 +4710,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) + & *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) + & 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) + & 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) + & 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 @@ -4800,9 +4812,9 @@ C fac_shield(i)=0.4d0 C fac_shield(j)=0.6d0 endif ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss C Diagnostics. Comment out or remove after debugging! c ees0p(num_conti,i)=0.5D0*fac3*ees0pij c ees0m(num_conti,i)=0.5D0*fac3*ees0mij @@ -4851,11 +4863,17 @@ cd fprimcont=0.0D0 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k) enddo gggp(1)=gggp(1)+ees0pijp*xj + & +ees0p(num_conti,i)/sss*rmij*xj*sssgrad gggp(2)=gggp(2)+ees0pijp*yj + & +ees0p(num_conti,i)/sss*rmij*yj*sssgrad gggp(3)=gggp(3)+ees0pijp*zj + & +ees0p(num_conti,i)/sss*rmij*zj*sssgrad gggm(1)=gggm(1)+ees0mijp*xj + & +ees0m(num_conti,i)/sss*rmij*xj*sssgrad gggm(2)=gggm(2)+ees0mijp*yj + & +ees0m(num_conti,i)/sss*rmij*yj*sssgrad gggm(3)=gggm(3)+ees0mijp*zj + & +ees0m(num_conti,i)/sss*rmij*zj*sssgrad C Derivatives due to the contact function gacont_hbr(1,num_conti,i)=fprimcont*xj gacont_hbr(2,num_conti,i)=fprimcont*yj @@ -4870,28 +4888,28 @@ cgrad ghalfm=0.5D0*gggm(k) gacontp_hb1(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontp_hb2(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontp_hb3(k,num_conti,i)=gggp(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb1(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb2(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb3(k,num_conti,i)=gggm(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss enddo C Diagnostics. Comment out or remove after debugging! @@ -4957,6 +4975,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 @@ -4994,7 +5014,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, @@ -5003,10 +5023,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 @@ -5061,14 +5081,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) @@ -5082,7 +5102,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) @@ -5091,7 +5111,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 @@ -5099,7 +5119,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) @@ -5107,8 +5127,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------------------------------------------------------------------------------- @@ -5237,7 +5266,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) @@ -5285,12 +5314,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 @@ -5320,7 +5343,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)) @@ -5329,7 +5352,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)) @@ -5341,7 +5364,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 @@ -5361,7 +5384,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 @@ -5380,7 +5403,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) @@ -5395,7 +5418,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) @@ -5410,7 +5433,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) @@ -5426,8 +5449,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----------------------------------------------------------------------------- @@ -5482,7 +5513,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CONTROL' dimension ggg(3) - integer xshift,yshift,zshift + double precision boxshift evdw2=0.0D0 evdw2_14=0.0d0 r0_scp=4.5d0 @@ -5491,7 +5522,10 @@ 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 +c do i=iatscp_s,iatscp_e + do ikont=g_listscp_start,g_listscp_end + i=newcontlistscpi(ikont) + j=newcontlistscpj(ikont) 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)) @@ -5522,18 +5556,13 @@ 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 + call to_box(xi,yi,zi) C xi=xi+xshift*boxxsize C yi=yi+yshift*boxysize C zi=zi+zshift*boxzsize - do iint=1,nscp_gr(i) +c do iint=1,nscp_gr(i) - do j=iscpstart(i,iint),iscpend(i,iint) +c do j=iscpstart(i,iint),iscpend(i,iint) if (itype(j).eq.ntyp1) cycle itypj=iabs(itype(j)) C Uncomment following three lines for SC-p interactions @@ -5544,67 +5573,10 @@ C Uncomment following three lines for Ca-p interactions 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 + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) C xj=xj-xi C yj=yj-yi C zj=zj-zi @@ -5626,39 +5598,13 @@ C ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac -cgrad if (j.lt.i) then -cd write (iout,*) 'ji' -cgrad do k=1,3 -cgrad ggg(k)=-ggg(k) -C Uncomment following line for SC-p interactions -c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k) -cgrad enddo -cgrad endif -cgrad do k=1,3 -cgrad gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k) -cgrad enddo -cgrad kstart=min0(i+1,j) -cgrad kend=max0(i-1,j-1) -cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend -cd write (iout,*) ggg(1),ggg(2),ggg(3) -cgrad do k=kstart,kend -cgrad do l=1,3 -cgrad gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l) -cgrad enddo -cgrad enddo do k=1,3 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - enddo +c enddo - enddo ! iint +c enddo ! iint enddo ! i C enddo !zshift C enddo !yshift @@ -5672,7 +5618,10 @@ C This subroutine calculates the excluded-volume interaction energy between C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C - implicit real*8 (a-h,o-z) + implicit none +#ifdef MPI + include 'mpif.h' +#endif include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -5684,8 +5633,18 @@ C include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.SPLITELE' - integer xshift,yshift,zshift - dimension ggg(3) + include 'COMMON.TIME1' + double precision ggg(3) + integer i,iint,j,k,iteli,itypj,subchap,ikont + double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, + & fac,e1,e2,rij + double precision evdw2,evdw2_14,evdwij + double precision sscale,sscagrad + double precision boxshift + external boxshift,to_box +c#ifdef TIMING_ENE +c double precision time01 +c#endif evdw2=0.0D0 evdw2_14=0.0d0 c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' @@ -5694,53 +5653,24 @@ 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 - if (energy_dec) write (iout,*) "escp:",r_cut,rlamb - do i=iatscp_s,iatscp_e + if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb +c do i=iatscp_s,iatscp_e + do ikont=g_listscp_start,g_listscp_end +c#ifdef TIMING_ENE +c time01=MPI_Wtime() +c#endif + i=newcontlistscpi(ikont) + j=newcontlistscpj(ikont) 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) +!DIR$ INLINE + call to_box(xi,yi,zi) +c do iint=1,nscp_gr(i) - do j=iscpstart(i,iint),iscpend(i,iint) +c do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions @@ -5751,76 +5681,28 @@ C Uncomment following three lines for Ca-p interactions 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 +!DIR$ INLINE + call to_box(xj,yj,zj) +c#ifdef TIMING_ENE +c time_escpsetup=time_escpsetup+MPI_Wtime()-time01 +c time01=MPI_Wtime() +c#endif +!DIR$ INLINE + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) c print *,xj,yj,zj,'polozenie j' +c#ifdef TIMING_ENE +c time_escpsetup=time_escpsetup+MPI_Wtime()-time01 +c time01=MPI_Wtime() +c#endif rrij=1.0D0/(xj*xj+yj*yj+zj*zj) c print *,rrij - sss=sscale(1.0d0/(dsqrt(rrij))) + sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int) 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))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij)),r_cut_int) fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) @@ -5831,8 +5713,9 @@ c if (sss.eq.0) print *,'czasem jest OK' endif evdwij=e1+e2 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), + if (energy_dec) write (iout,'(a6,2i5,3f7.3,2i3,3e11.3)') + & 'evdw2',i,j,1.0d0/dsqrt(rrij),sss, + & evdwij,iteli,itypj,fac,aad(itypj,iteli), & bad(itypj,iteli) C C Calculate contributions to the gradient in the virtual-bond and SC vectors. @@ -5873,10 +5756,13 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo +c#ifdef TIMING_ENE +c time_escpcalc=time_escpcalc+MPI_Wtime()-time01 +c#endif c endif !endif for sscale cutoff - enddo ! j +c enddo ! j - enddo ! iint +c enddo ! iint enddo ! i c enddo !zshift c enddo !yshift @@ -6058,7 +5944,8 @@ C 15/02/13 CC dynamic SSbond - additional check if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) - ehpb=ehpb+2*eij +c ehpb=ehpb+2*eij + ehpb=ehpb+eij endif cd write (iout,*) "eij",eij cd & ' waga=',waga,' fac=',fac @@ -7224,7 +7111,8 @@ c & sumene4, c & dscp1,dscp2,sumene c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene -c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i) + if (energy_dec) write (2,*) "i",i," itype",itype(i)," it",it, + & " escloc",sumene,escloc,it,itype(i) c & ,zz,xx,yy c#define DEBUG #ifdef DEBUG @@ -8083,16 +7971,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 @@ -12170,6 +12062,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 @@ -12224,6 +12117,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) @@ -12234,6 +12129,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 @@ -13536,3 +13433,111 @@ C----------------------------------------------------------------------- endif return end +c------------------------------------------------------------------------ + double precision function boxshift(x,boxsize) + implicit none + double precision x,boxsize + double precision xtemp + xtemp=dmod(x,boxsize) + if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then + boxshift=xtemp-boxsize + else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then + boxshift=xtemp+boxsize + else + boxshift=xtemp + endif + return + end +c-------------------------------------------------------------------------- + subroutine closest_img(xi,yi,zi,xj,yj,zj) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + integer xshift,yshift,zshift,subchap + double precision dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp + xj_safe=xj + yj_safe=yj + zj_safe=zj + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + 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 + return + end +c-------------------------------------------------------------------------- + subroutine to_box(xi,yi,zi) + implicit none + include 'DIMENSIONS' + include 'COMMON.CHAIN' + double precision xi,yi,zi + xi=dmod(xi,boxxsize) + if (xi.lt.0.0d0) xi=xi+boxxsize + yi=dmod(yi,boxysize) + if (yi.lt.0.0d0) yi=yi+boxysize + zi=dmod(zi,boxzsize) + if (zi.lt.0.0d0) zi=zi+boxzsize + return + end +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 +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 +#ifdef DEBUG + write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi +#endif + return + end