C Calculate the virtual-bond-angle energy.
C
if (wang.gt.0d0) then
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
else
ebe=0
endif
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"
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
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)
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
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
#ifdef DEBUG
write (iout,*) "gloc after reduce"
do i=1,4*nres
& 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)'/
& '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)')
& 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)'/
& '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)')
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
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))
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
& 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--------------------------------------------------------------------------
& 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
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
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
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
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp((itype(i)))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
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)
& 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
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,
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
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)
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
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
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,
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
& auxvec1(2),auxmat1(2,2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C C
+C C
C Parallel Antiparallel C
C C
C o o 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
projects / unres.git / blobdiff