X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=ec0aef92c147b72de51eb82ba118958ee2a6d49f;hb=42b30ef61fe3a8077aee8292bd43e4de4e6b8522;hp=5b7cc93a64b9f861dcc26121d8fa748eb218eaf9;hpb=e3a9e82f4769b262c53a0b3d0aa297880600d179;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 5b7cc93..ec0aef9 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -190,7 +190,7 @@ C C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then - call ebend(ebe) + call ebend(ebe,ethetacnstr) else ebe=0 endif @@ -301,6 +301,7 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor +C energia(22)= c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" @@ -716,7 +717,7 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -725,7 +726,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG do i=1,nres do j=1,3 gloc_scbuf(j,i)=gloc_sc(j,i,icg) @@ -745,7 +746,7 @@ c enddo call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -754,7 +755,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -977,7 +978,8 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, - & edihcnstr,ebr*nss, + & edihcnstr, + & ethetacnstr,ebr*nss, & Uconst,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ @@ -1000,6 +1002,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & 'ETOT= ',1pE16.6,' (total)') @@ -1010,6 +1013,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, + & ethetacnstr, & ebr*nss,Uconst,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ @@ -1031,6 +1035,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & 'ETOT= ',1pE16.6,' (total)') @@ -1420,6 +1425,9 @@ C include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' logical lprn + +c write(iout,*) "Jestem w egb(evdw)" + evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -1447,10 +1455,36 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + +c write(iout,*) "PRZED ZWYKLE", evdwij call dyn_ssbond_ene(i,j,evdwij) +c write(iout,*) "PO ZWYKLE", evdwij + 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 search over all next residues + 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) +c if(evdwij_przed_tri.ne.evdwij) then +c write (iout,*) "TRI:", evdwij, evdwij_przed_tri +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',i,j,evdwij,'tss' + endif!dyn_ss_mask(k) + enddo! k ELSE ind=ind+1 itypj=iabs(itype(j)) @@ -4055,8 +4089,13 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 + do i=1,3 + ggg(i)=0.0d0 + enddo +C write (iout,*) ,"link_end",link_end,constr_dist cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr cd write(iout,*)'link_start=',link_start,' link_end=',link_end if (link_end.eq.0) return @@ -4085,47 +4124,104 @@ C 15/02/13 CC dynamic SSbond - additional check & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij -cd write (iout,*) "eij",eij endif +cd write (iout,*) "eij",eij +cd & ' waga=',waga,' fac=',fac + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, + & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo else C Calculate the distance between the two points and its difference from the C target distance. - dd=dist(ii,jj) - rdis=dd-dhpb(i) + dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, + & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. - waga=forcon(i) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C - fac=waga*rdis/dd -cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, -cd & ' waga=',waga,' fac=',fac - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo + fac=waga*rdis/dd + endif + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). - if (iii.lt.ii) then + if (iii.lt.ii) then do j=1,3 ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo - endif + endif cgrad do j=iii,jjj-1 cgrad do k=1,3 cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k) cgrad enddo cgrad enddo - do k=1,3 - ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) - ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) - enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -4429,6 +4525,31 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai + do i=ithetaconstr_start,ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itori) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif + enddo if (energy_dec) write (iout,'(a6,i5,0pf7.3)') & 'ebend',i,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 @@ -4574,7 +4695,9 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).eq.ntyp1)) cycle +C print *,i,theta(i) if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -4586,7 +4709,9 @@ C coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then +C print *,ethetai + + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4601,13 +4726,13 @@ C propagation of chirality for glycine type enddo else phii=0.0d0 - ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4622,7 +4747,7 @@ C propagation of chirality for glycine type enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -4672,6 +4797,7 @@ C propagation of chirality for glycine type enddo write(iout,*) "ethetai",ethetai endif +C print *,ethetai do m=1,ntheterm2 do k=1,nsingle aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k) @@ -4692,10 +4818,16 @@ C propagation of chirality for glycine type & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet", & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet", & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) enddo enddo +C print *,"cosph1", (cosph1(k), k=1,nsingle) +C print *,"cosph2", (cosph2(k), k=1,nsingle) +C print *,"sinph1", (sinph1(k), k=1,nsingle) +C print *,"sinph2", (sinph2(k), k=1,nsingle) if (lprn) & write(iout,*) "ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) do m=1,ntheterm3 do k=2,ndouble do l=1,k-1 @@ -4730,7 +4862,35 @@ C propagation of chirality for glycine type enddo enddo 10 continue +C now we have the theta_constrains + do i=ithetaconstr_start,ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itori) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif + enddo + c lprn1=.true. +C print *,ethetai if (lprn1) & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & i,theta(i)*rad2deg,phii*rad2deg, @@ -4739,7 +4899,7 @@ c lprn1=.false. etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo return end @@ -5560,12 +5720,12 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) @@ -5664,18 +5824,21 @@ c do i=1,ndih_constr difi=pinorm(phii-phi0(i)) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else difi=0.0 endif -cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, -cd & rad2deg*phi0(i), rad2deg*drange(i), -cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc", + & i,itori,rad2deg*phii, + & rad2deg*phi0(i), rad2deg*drange(i), + & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) + endif enddo cd write (iout,*) 'edihcnstr',edihcnstr return @@ -7733,9 +7896,9 @@ cd ghalf=0.0d0 cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k) cgrad ghalf=0.5d0*ggg2(ll) cd ghalf=0.0d0 - gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2) + gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2) gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2) - gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2) + gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2) gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2) gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl @@ -8040,12 +8203,12 @@ C C C o o C C \ /l\ /j\ / C C \ / \ / \ / C -C o| o | | o |o C +C o| o | | o |o C C \ j|/k\| \ |/k\|l C C \ / \ \ / \ C C o o C -C i i C -C C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l C AL 7/4/01 s1 would occur in the sixth-order moment, @@ -8216,10 +8379,10 @@ c---------------------------------------------------------------------------- double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C C +C C C Parallel Antiparallel C C C -C o o C +C o o C C /l\ / \ /j\ C C / \ / \ / \ C C /| o |o o| o |\ C @@ -8333,7 +8496,7 @@ c---------------------------------------------------------------------------- & auxvec1(2),auxmat1(2,2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C C +C C C Parallel Antiparallel C C C C o o C @@ -8341,10 +8504,10 @@ C /l\ / \ /j\ C C / \ / \ / \ C C /| o |o o| o |\ C C \ j|/k\| \ |/k\|l C -C \ / \ \ / \ C +C \ / \ \ / \ C C o \ o \ C C i i C -C C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective