X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=ddc9833cd43ff65c32bf3957f2752600f8acda9e;hb=bc23440fbe68672d430f71f22f46b11265f003db;hp=a39509acef09f8a3e549eb4816344255fa2cd338;hpb=3d3780d89d2885ee8700324ef832ef07ee35c45b;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index a39509a..ddc9833 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -400,12 +400,14 @@ if (nfgtasks.gt.1) then call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR) endif + if (nres_molec(1).gt.0) then if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list ! write (iout,*) "after make_SCp_inter_list" if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list ! write (iout,*) "after make_SCSC_inter_list" if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list + endif ! write (iout,*) "after make_pp_inter_list" ! print *,'Processor',myrank,' calling etotal ipot=',ipot @@ -1930,7 +1932,7 @@ dGCLdOM1=0.0d0 dPOLdOM1=0.0d0 ! write (iout,*) "RWA", g_listscsc_start,g_listscsc_end,i,j - + if (nres_molec(1).eq.0) return do icont=g_listscsc_start,g_listscsc_end i=newcontlisti(icont) j=newcontlistj(icont) @@ -1966,7 +1968,7 @@ 'evdw',i,j,evdwij,' ss' ! if (energy_dec) write (iout,*) & ! 'evdw',i,j,evdwij,' ss' - do k=j+1,iend(i,iint) + do k=j+1,nres !C search over all next residues if (dyn_ss_mask(k)) then !C check if they are cysteins @@ -3463,6 +3465,7 @@ eel_loc=0.0d0 eello_turn3=0.0d0 eello_turn4=0.0d0 + if (nres_molec(1).eq.0) return ! if (icheckgrad.eq.1) then @@ -5076,14 +5079,19 @@ #ifdef NEWCORR gloc(nphi+i,icg)=gloc(nphi+i,icg)& -(gs13+gsE13+gsEE1)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)& -(gs23+gs21+gsEE2)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j)& + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)& -(gs32+gsE31+gsEE3)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j)& + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + !c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- !c & gs2 @@ -5409,6 +5417,7 @@ !d print '(a)','Enter ESCP' !d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e ! do i=iatscp_s,iatscp_e + if (nres_molec(1).eq.0) return do icont=g_listscp_start,g_listscp_end i=newcontlistscpi(icont) j=newcontlistscpj(icont) @@ -5570,7 +5579,7 @@ iabs(itype(jjj,1)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij -!d write (iout,*) "eij",eij +! write (iout,*) "eij",eij,iii,jjj endif else if (ii.gt.nres .and. jj.gt.nres) then !c Restraints from contact prediction @@ -5709,10 +5718,13 @@ itypj=iabs(itype(j,1)) ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(nres+j) - 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) + 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) @@ -5737,9 +5749,9 @@ eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) & +akct*deltad*deltat12 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr -! write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, -! & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, -! & " deltat12",deltat12," eij",eij +! write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, & +! " akct",akct," deltad",deltad," deltat",deltat1,deltat2, & +! " deltat12",deltat12," eij",eij ed=2*akcm*deltad+akct*deltat12 pom1=akct*deltad pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi @@ -12436,7 +12448,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! enddo do k=1,3 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)) -!C print *,'gg',k,gg(k) +! print *,'gg',k,gg(k) enddo ! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut ! write (iout,*) "gg",(gg(k),k=1,3) @@ -13194,6 +13206,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! do j=1,3 grad_s(j,i)=gcart(j,i) grad_s(j+3,i)=gxcart(j,i) + write(iout,*) "before movement analytical gradient" + do i=1,nres + write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3) + enddo + enddo enddo else @@ -13421,12 +13439,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) -! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i) - -! if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo + write(iout,*) "before movement analytical gradient" + do i=1,nres + write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3) + enddo + else !- split gradient check call zerograd @@ -14095,6 +14116,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi call to_box(xj,yj,zj) + 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 @@ -17401,13 +17426,13 @@ chip1=chip(itypi) #endif !#define DEBUG !el write (iout,*) "After sum_gradient" -#ifdef DEBUG - write (iout,*) "After sum_gradient" - do i=1,nres-1 - write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) - write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) - enddo -#endif +!#ifdef DEBUG +! write (iout,*) "After sum_gradient" +! do i=1,nres-1 +! write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) +! write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) +! enddo +!#endif !#undef DEBUG ! If performing constraint dynamics, add the gradients of the constraint energy if(usampl.and.totT.gt.eq_time) then @@ -17446,8 +17471,8 @@ chip1=chip(itypi) ! if (i.le.2) print *,"gcart_one",gcart(j,i),gradc(j,i,icg) enddo #ifdef DEBUG - write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),& - (gxcart(j,i),j=1,3),gloc(i,icg) + write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3),gloc(i,icg),(gloc_sc(j,i,icg),j=1,3) #endif enddo #ifdef TIMING @@ -17778,10 +17803,12 @@ chip1=chip(itypi) do j=1,3 dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/& vbld(i-1) - if (itype(i-1,1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint + if (((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))) & + dtheta(j,1,i)=-dcostheta(j,1,i)/sint dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/& vbld(i) - if (itype(i-1,1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint + if ((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))& + dtheta(j,2,i)=-dcostheta(j,2,i)/sint enddo enddo #if defined(MPI) && defined(PARINTDER) @@ -17838,11 +17865,23 @@ chip1=chip(itypi) cost1=dcos(theta(i-1)) cosg=dcos(phi(i)) scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. & phi(i).gt.pi34.and.phi(i).le.pi.or. & @@ -17851,8 +17890,16 @@ chip1=chip(itypi) call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2) call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3) do j=1,3 + if (sint.ne.0.0d0) then ctgt=cost/sint + else + ctgt=0.0d0 + endif + if (sint1.ne.0.0d0) then ctgt1=cost1/sint1 + else + ctgt1=0.0d0 + endif cosg_inv=1.0d0/cosg if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) & @@ -17867,6 +17914,10 @@ chip1=chip(itypi) ! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1) dphi(j,3,i)=cosg_inv*dsinphi(j,3,i) endif +! write(iout,*) "just after,close to pi",dphi(j,3,i),& +! sing*(ctgt1*dtheta(j,2,i-1)),ctgt*dtheta(j,1,i), & +! (fac0*vp2(j)+sing*dc_norm(j,i-2)),vbld_inv(i-1) + ! Bug fixed 3/24/05 (AL) enddo ! Obtaining the gamma derivatives from cosine derivative @@ -17921,12 +17972,25 @@ chip1=chip(itypi) ! write(iout,*) dc_norm2(j,i-2+nres),"dcnorm" enddo scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1)) + ! write(iout,*) "faki",fac0,fac1,fac2,fac3,fac + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) - ! write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4 + else + fac4=0.0d0 + endif + ! Obtaining the gamma derivatives from sine derivative if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or. & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or. & @@ -17994,11 +18058,23 @@ chip1=chip(itypi) ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) ! enddo scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or. & tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or. & @@ -18069,11 +18145,23 @@ chip1=chip(itypi) ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) enddo scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or. & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or. & @@ -19731,14 +19819,12 @@ chip1=chip(itypi) if (idssb(i).eq.newihpb(j) .and. & jdssb(i).eq.newjhpb(j)) found=.true. enddo -#ifndef CLUST -#ifndef WHAM +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) ! write(iout,*) "found",found,i,j if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i) #endif -#endif enddo do i=1,newnss @@ -19748,21 +19834,22 @@ chip1=chip(itypi) if (newihpb(i).eq.idssb(j) .and. & newjhpb(i).eq.jdssb(j)) found=.true. enddo -#ifndef CLUST -#ifndef WHAM +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) ! write(iout,*) "found",found,i,j if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i) #endif -#endif enddo - +!#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) nss=newnss do i=1,nss idssb(i)=newihpb(i) jdssb(i)=newjhpb(i) enddo +!#else +! nss=0 +!#endif return end subroutine dyn_set_nss @@ -21242,10 +21329,10 @@ chip1=chip(itypi) !(3,3,2,maxres) ! allocateion of lists JPRDLA allocate(newcontlistppi(300*nres)) - allocate(newcontlistscpi(300*nres)) + allocate(newcontlistscpi(350*nres)) allocate(newcontlisti(300*nres)) allocate(newcontlistppj(300*nres)) - allocate(newcontlistscpj(300*nres)) + allocate(newcontlistscpj(350*nres)) allocate(newcontlistj(300*nres)) return @@ -22969,6 +23056,10 @@ chip1=chip(itypi) real(kind=8) :: facd4, adler, Fgb, facd3 integer troll,jj,istate real (kind=8) :: dcosom1(3),dcosom2(3) + real(kind=8) ::locbox(3) + locbox(1)=boxxsize + locbox(2)=boxysize + locbox(3)=boxzsize evdw=0.0D0 if (nres_molec(5).eq.0) return @@ -23011,6 +23102,8 @@ chip1=chip(itypi) zj=c(3,j) call to_box(xj,yj,zj) +! write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,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 @@ -23019,6 +23112,7 @@ chip1=chip(itypi) xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) +! write(iout,*) "xj,yj,zj", xj,yj,zj,boxxsize ! dxj = dc_norm( 1, nres+j ) ! dyj = dc_norm( 2, nres+j ) @@ -23051,6 +23145,11 @@ chip1=chip(itypi) b3cav = alphasurcat(3,itypi,itypj) b4cav = alphasurcat(4,itypi,itypj) +! b1cav=0.0d0 +! b2cav=0.0d0 +! b3cav=0.0d0 +! b4cav=0.0d0 + ! used to determine whether we want to do quadrupole calculations eps_in = epsintabcat(itypi,itypj) if (eps_in.eq.0.0) eps_in=1.0 @@ -23062,11 +23161,13 @@ chip1=chip(itypi) ctail(k,1)=c(k,i+nres) ctail(k,2)=c(k,j) END DO + call to_box(ctail(1,1),ctail(2,1),ctail(3,1)) + call to_box(ctail(1,2),ctail(2,2),ctail(3,2)) !c! tail distances will be themselves usefull elswhere !c1 (in Gcav, for example) - Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 ) - Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 ) - Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 ) + do k=1,3 + Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k)) + enddo Rtail = dsqrt( & (Rtail_distance(1)*Rtail_distance(1)) & + (Rtail_distance(2)*Rtail_distance(2)) & @@ -23081,10 +23182,15 @@ chip1=chip(itypi) ! see unres publications for very informative images chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres) chead(k,2) = c(k, j) + enddo + call to_box(chead(1,1),chead(2,1),chead(3,1)) + call to_box(chead(1,2),chead(2,2),chead(3,2)) +! write(iout,*) "TEST",chead(1,1),chead(2,1),chead(3,1),dc_norm(k, i+nres),d1 ! distance ! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) -! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) - Rhead_distance(k) = chead(k,2) - chead(k,1) +! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) + do k=1,3 + Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k)) END DO ! pitagoras (root of sum of squares) Rhead = dsqrt( & @@ -23106,6 +23212,7 @@ chip1=chip(itypi) dPOLdOM1 = 0.0d0 dPOLdOM2 = 0.0d0 Fcav = 0.0d0 + Fisocav=0.0d0 dFdR = 0.0d0 dCAVdOM1 = 0.0d0 dCAVdOM2 = 0.0d0 @@ -23131,6 +23238,16 @@ chip1=chip(itypi) rij_shift = Rtail - sig + sig0ij IF (rij_shift.le.0.0D0) THEN evdw = 1.0D20 + if (evdw.gt.1.0d6) then + write (*,'(2(1x,a3,i3),7f7.2)') & + restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& + 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq + write(*,*) facsig,faceps1_inv,om1,chiom1,chi1 + write(*,*) "ANISO?!",chi1 +!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw + endif + RETURN END IF sigder = -sig * sigsq @@ -23164,7 +23281,7 @@ chip1=chip(itypi) gg(1) = fac gg(2) = fac gg(3) = fac - +! print *,"GG(1),distance grad",gg(1) fac = chis1 * sqom1 + chis2 * sqom2 & - 2.0d0 * chis12 * om1 * om2 * om12 pom = 1.0d0 - chis1 * chis2 * sqom12 @@ -23203,12 +23320,12 @@ chip1=chip(itypi) erdxi = scalar( ertail(1), dC_norm(1,i+nres) ) erdxj = scalar( ertail(1), dC_norm(1,j) ) facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) - facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres) + facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j) DO k = 1, 3 pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres)) gradpepcatx(k,i) = gradpepcatx(k,i) & - (( dFdR + gg(k) ) * pom) - pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres)) + pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j)) ! gvdwx(k,j) = gvdwx(k,j) & ! + (( dFdR + gg(k) ) * pom) gradpepcat(k,i) = gradpepcat(k,i) & @@ -23218,7 +23335,10 @@ chip1=chip(itypi) gg(k) = 0.0d0 ENDDO !c! Compute head-head and head-tail energies for each state +!! if (.false.) then ! turn off electrostatic + if (itype(j,5).gt.0) then ! the normal cation case isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj) +! print *,i,itype(i,1),isel IF (isel.eq.0) THEN !c! No charges - do nothing eheadtail = 0.0d0 @@ -23229,10 +23349,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL enq_cat(epol) eheadtail = epol @@ -23244,10 +23360,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif ! write(iout,*) "KURWA0",d1 CALL edq_cat(ecl, elj, epol) @@ -23261,10 +23373,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL eqq_cat(Ecl,Egb,Epol,Fisocav,Elj) eheadtail = ECL + Egb + Epol + Fisocav + Elj @@ -23285,16 +23393,28 @@ chip1=chip(itypi) ! ! CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav + else + write(iout,*) "not yet implemented",j,itype(j,5) + endif +!! endif ! turn off electrostatic evdw = evdw + Fcav + eheadtail +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& Equad,evdwij+Fcav+eheadtail,evdw ! evdw = evdw + Fcav + eheadtail - +! print *,"before sc_grad_cat", i,j, gradpepcat(1,j) ! iF (nstate(itypi,itypj).eq.1) THEN CALL sc_grad_cat +! print *,"after sc_grad_cat", i,j, gradpepcat(1,j) + ! END IF !c!------------------------------------------------------------------- !c! NAPISY KONCOWE @@ -23305,6 +23425,7 @@ chip1=chip(itypi) ! print *,"EVDW KURW",evdw,nres !!! return 17 continue +! go to 23 do i=ibond_start,ibond_end ! print *,"I am in EVDW",i @@ -23333,6 +23454,10 @@ chip1=chip(itypi) yj=c(2,j) zj=c(3,j) call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 dxj = 0.0d0! dc_norm( 1, nres+j ) @@ -23377,11 +23502,16 @@ chip1=chip(itypi) ctail(k,1)=(c(k,i)+c(k,i+1))/2.0 ctail(k,2)=c(k,j) END DO + call to_box(ctail(1,1),ctail(2,1),ctail(3,1)) + call to_box(ctail(1,2),ctail(2,2),ctail(3,2)) +!c! tail distances will be themselves usefull elswhere +!c1 (in Gcav, for example) + do k=1,3 + Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k)) + enddo + !c! tail distances will be themselves usefull elswhere !c1 (in Gcav, for example) - Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 ) - Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 ) - Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 ) Rtail = dsqrt( & (Rtail_distance(1)*Rtail_distance(1)) & + (Rtail_distance(2)*Rtail_distance(2)) & @@ -23398,11 +23528,20 @@ chip1=chip(itypi) ! see unres publications for very informative images chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i) chead(k,2) = c(k, j) + ENDDO ! distance ! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) ! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) - Rhead_distance(k) = chead(k,2) - chead(k,1) + call to_box(chead(1,1),chead(2,1),chead(3,1)) + call to_box(chead(1,2),chead(2,2),chead(3,2)) + +! distance +! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) +! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) + do k=1,3 + Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k)) END DO + ! pitagoras (root of sum of squares) Rhead = dsqrt( & (Rhead_distance(1)*Rhead_distance(1)) & @@ -23448,6 +23587,13 @@ chip1=chip(itypi) rij_shift = Rtail - sig + sig0ij IF (rij_shift.le.0.0D0) THEN evdw = 1.0D20 +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),6f6.2)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij +!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif RETURN END IF sigder = -sig * sigsq @@ -23539,24 +23685,26 @@ chip1=chip(itypi) + (( dFdR + gg(k) ) * ertail(k)) gg(k) = 0.0d0 ENDDO + if (itype(j,5).gt.0) then !c! Compute head-head and head-tail energies for each state isel = 3 !c! Dipole-charge interactions - if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then - Qi=Qi*2 - Qij=Qij*2 - endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL edq_cat_pep(ecl, elj, epol) eheadtail = ECL + elj + epol ! print *,"i,",i,eheadtail ! eheadtail = 0.0d0 - + else +!HERE WATER and other types of molecules solvents will be added + write(iout,*) "not yet implemented" +! CALL edd_cat_pep + endif evdw = evdw + Fcav + eheadtail - +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& @@ -23573,7 +23721,8 @@ chip1=chip(itypi) !c write (iout,*) "Number of loop steps in EGB:",ind !c energy_dec=.false. ! print *,"EVDW KURW",evdw,nres - + 23 continue +! print *,"before leave sc_grad_cat", i,j, gradpepcat(1,nres-1) return end subroutine ecats_prot_amber @@ -24128,9 +24277,13 @@ chip1=chip(itypi) dEvan2Cm,cm1,cm,vcat,vsug,v1,v2,dx,vcm,dEdipCm,dEdipCalp, & dEvan1Calp,dEvan2Cat,dEvan2Calp,dEtotalCat,dEdipCat,dEvan1Cat,dcosdcat, & dcosdcalp,dcosdcm,dEgbdCat,dEgbdCalp,dEgbdCm,dEcavdCat,dEcavdCalp, & - dEcavdCm + dEcavdCm,boxik real(kind=8),dimension(14) :: vcatnuclprm ecation_nucl=0.0d0 + boxik(1)=boxxsize + boxik(2)=boxysize + boxik(3)=boxzsize + if (nres_molec(5).eq.0) return itmp=0 do i=1,4 @@ -24151,6 +24304,7 @@ chip1=chip(itypi) yj=c(2,j) zj=c(3,j) call to_box(xj,yj,zj) +! write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,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 @@ -24159,7 +24313,7 @@ chip1=chip(itypi) xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) - +! write(iout,*) 'after shift', xj,yj,zj dist_init=xj**2+yj**2+zj**2 itypi=itype(i,2) @@ -24172,12 +24326,16 @@ chip1=chip(itypi) vsug(k)=c(k,i) vcat(k)=c(k,j) enddo + call to_box(vcm(1),vcm(2),vcm(3)) + call to_box(vsug(1),vsug(2),vsug(3)) + call to_box(vcat(1),vcat(2),vcat(3)) do k=1,3 - dx(k) = vcat(k)-vcm(k) - enddo - do k=1,3 +! dx(k) = vcat(k)-vcm(k) +! enddo + dx(k)=boxshift(vcat(k)-vcm(k),boxik(k)) +! do k=1,3 v1(k)=dc(k,i+nres) - v2(k)=(vcat(k)-vsug(k)) + v2(k)=boxshift(vcat(k)-vsug(k),boxik(k)) enddo v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2) v1dpdx = v1(1)*dx(1)+v1(2)*dx(2)+v1(3)*dx(3) @@ -25938,28 +26096,7 @@ chip1=chip(itypi) zi=c(3,nres+i) call to_box(xi,yi,zi) call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) - 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 +! endif ! print *, sslipi,ssgradlipi dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) @@ -26016,7 +26153,7 @@ chip1=chip(itypi) zj=c(3,j+nres) call to_box(xj,yj,zj) call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) - write(iout,*) "KRUWA", i,j +! write(iout,*) "KRUWA", i,j 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 & @@ -26594,7 +26731,7 @@ chip1=chip(itypi) ! integer :: k !c! Epol and Gpol analytical parameters alphapol1 = alphapolcat(itypi,itypj) - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat2(itypj,itypi) !c! Fisocav and Gisocav analytical parameters al1 = alphisocat(1,itypi,itypj) al2 = alphisocat(2,itypi,itypj) @@ -27247,7 +27384,7 @@ chip1=chip(itypi) use calc_data use comm_momo double precision facd3, adler,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) !c! R2 - distance between head of jth side chain and tail of ith sidechain R2 = 0.0d0 DO k = 1, 3 @@ -27545,7 +27682,7 @@ chip1=chip(itypi) use calc_data double precision facd3, adler,ecl,elj,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) w1 = wqdipcat(1,itypi,itypj) w2 = wqdipcat(2,itypi,itypj) pis = sig0headcat(itypi,itypj) @@ -27664,7 +27801,7 @@ chip1=chip(itypi) use calc_data double precision facd3, adler,ecl,elj,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) w1 = wqdipcat(1,itypi,itypj) w2 = wqdipcat(2,itypi,itypj) pis = sig0headcat(itypi,itypj) @@ -28010,9 +28147,9 @@ chip1=chip(itypi) alf1 = 0.0d0 alf2 = 0.0d0 alf12 = 0.0d0 - dxj = dc_norm( 1, nres+j ) - dyj = dc_norm( 2, nres+j ) - dzj = dc_norm( 3, nres+j ) + dxj = 0.0d0 !dc_norm( 1, nres+j ) + dyj = 0.0d0 !dc_norm( 2, nres+j ) + dzj = 0.0d0 !dc_norm( 3, nres+j ) !c! distance from center of chain(?) to polar/charged head d1 = dheadcat(1, 1, itypi, itypj) d2 = dheadcat(2, 1, itypi, itypj) @@ -28262,8 +28399,10 @@ chip1=chip(itypi) zi=c(3,nres+i) call to_box(xi,yi,zi) do iint=1,nint_gr(i) +! print *,"is it wrong", iint,i do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j,1)) + if (energy_dec) write(iout,*) "LISTA ZAKRES",istart(i,iint),iend(i,iint),iatsc_s,iatsc_e if (itypj.eq.ntyp1) cycle xj=c(1,nres+j) yj=c(2,nres+j) @@ -28350,7 +28489,7 @@ chip1=chip(itypi) include 'mpif.h' real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp real*8 :: dist_init, dist_temp,r_buff_list - integer:: contlistscpi(250*nres),contlistscpj(250*nres) + integer:: contlistscpi(350*nres),contlistscpj(350*nres) ! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres) integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr