C Calculate the virtual-bond-angle energy.
C
if (wang.gt.0d0) then
- call ebend(ebe)
- else
+ call ebend(ebe)
+ else
ebe=0
endif
c print *,"Processor",myrank," computed UB"
else
etors=0
edihcnstr=0
- endif
+ endif
if (constr_homology.ge.1) then
call e_modeller(ehomology_constr)
C 6/23/01 Calculate double-torsional energy
C
if (wtor_d.gt.0) then
- call etor_d(etors_d)
+ call etor_d(etors_d)
else
- etors_d=0
+ etors_d=0
endif
c print *,"Processor",myrank," computed Utord"
C
c write (iout,*) "From EsaxsC: Esaxs_constr",Esaxs_constr
else
Esaxs_constr = 0.0d0
- endif
+ endif
c print *,"Processor",myrank," computed Ucorr"
C
C If performing constraint dynamics, call the constraint energy
C after the equilibration time
if(usampl.and.totT.gt.eq_time) then
call EconstrQ
- call Econstr_back
+ call Econstr_back
else
Uconst=0.0d0
Uconst_back=0.0d0
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
#ifdef DEBUG
write (iout,*) "sum_gradient gsaxsc, gsaxsx"
do i=0,nres
- write (iout,'(i3,3e15.5,5x,3e15.5)')
+ write (iout,'(i3,3e15.5,5x,3e15.5)')
& i,(gsaxsc(j,i),j=1,3),(gsaxsx(j,i),j=1,3)
enddo
call flush(iout)
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& 0.5d0*(wscp*gvdwc_scpp(j,i)+
- & welec*gelc_long(j,i) +
+ & welec*gelc_long(j,i)+
& wel_loc*gel_loc_long(j,i)+
& wcorr*gcorr_long(j,i)+
& wcorr5*gradcorr5_long(j,i)+
include 'COMMON.SBRIDGE'
logical lprn
integer xshift,yshift,zshift
+
evdw=0.0D0
ccccc energy_dec=.false.
C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
c write(iout,*) 'b1=',b1(1,i-2)
c write (iout,*) 'theta=', theta(i-1)
- enddo
+ enddo
#else
if (i.gt. nnt+2 .and. i.lt.nct+2) then
iti = itortyp(itype(i-2))
b1(2,i-2)=b(5,iti)
b2(1,i-2)=b(2,iti)
b2(2,i-2)=b(4,iti)
- b1tilde(1,i-2)=b1(1,i-2)
- b1tilde(2,i-2)=-b1(2,i-2)
- b2tilde(1,i-2)=b2(1,i-2)
- b2tilde(2,i-2)=-b2(2,i-2)
+ b1tilde(1,i-2)= b1(1,i-2)
+ b1tilde(2,i-2)=-b1(2,i-2)
+ b2tilde(1,i-2)= b2(1,i-2)
+ b2tilde(2,i-2)=-b2(2,i-2)
EE(1,2,i-2)=eeold(1,2,iti)
EE(2,1,i-2)=eeold(2,1,iti)
EE(2,2,i-2)=eeold(2,2,iti)
cd write (iout,*) 'b1',b1(:,iti)
cd write (iout,*) 'b2',b2(:,iti)
cd write (iout,*) "phi(",i,")=",phi(i)," sin1",sin1," cos1",cos1
-cd write (iout,*) 'Ug',Ug(:,:,i-2)
+cd write (iout,*) 'Ug',Ug(:,:,i-2)
c if (i .gt. iatel_s+2) then
if (i .gt. nnt+2) then
call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2))
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
+
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
- dimension ggg(3)
+ include 'COMMON.CONTROL'
+ dimension ggg(3),ggg_peak(3,20)
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
+c write(iout,*)'link_start=',link_start,' link_end=',link_end,
+c & " constr_dist",constr_dist," link_start_peak",link_start_peak,
+c & " link_end_peak",link_end_peak
+ if (link_end.eq.0.and.link_end_peak.eq.0) return
+ if (link_end_peak.ne.0) then
+ do i=link_start_peak,link_end_peak
+ ehpb_peak=0.0d0
+c print *,"i",i," link_end_peak",link_end_peak," ipeak",
+c & ipeak(1,i),ipeak(2,i)
+ do ip=ipeak(1,i),ipeak(2,i)
+ ii=ihpb_peak(ip)
+ jj=jhpb_peak(ip)
+ dd=dist(ii,jj)
+ iip=ip-ipeak(1,i)+1
+C iii and jjj point to the residues for which the distance is assigned.
+ if (ii.gt.nres) then
+ iii=ii-nres
+ jjj=jj-nres
+ else
+ iii=ii
+ jjj=jj
+ endif
+ aux=rlornmr1(dd,dhpb_peak(ip),dhpb1_peak(ip),forcon_peak(ip))
+ aux=dexp(-scal_peak*aux)
+ ehpb_peak=ehpb_peak+aux
+ fac=rlornmr1prim(dd,dhpb_peak(ip),dhpb1_peak(ip),
+ & forcon_peak(ip))*aux/dd
+ do j=1,3
+ ggg_peak(j,iip)=fac*(c(j,jj)-c(j,ii))
+ enddo
+ if (energy_dec) write (iout,'(a6,3i5,6f10.3,i5)')
+ & "edisL",i,ii,jj,dd,dhpb_peak(ip),dhpb1_peak(ip),
+ & forcon_peak(ip),fordepth_peak(ip),ehpb_peak
+ enddo
+c write (iout,*) "ehpb_peak",ehpb_peak," scal_peak",scal_peak
+ ehpb=ehpb-fordepth_peak(ipeak(1,i))*dlog(ehpb_peak)/scal_peak
+ do ip=ipeak(1,i),ipeak(2,i)
+ iip=ip-ipeak(1,i)+1
+ do j=1,3
+ ggg(j)=ggg_peak(j,iip)/ehpb_peak
+ enddo
+ ii=ihpb_peak(ip)
+ jj=jhpb_peak(ip)
+C iii and jjj point to the residues for which the distance is assigned.
+ if (ii.gt.nres) then
+ iii=ii-nres
+ jjj=jj-nres
+ else
+ iii=ii
+ jjj=jj
+ endif
+ 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
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
+ enddo
+ enddo
+ endif
do i=link_start,link_end
C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
C CA-CA distance used in regularization of structure.
C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
if (.not.dyn_ss .and. i.le.nss) then
C 15/02/13 CC dynamic SSbond - additional check
- if (ii.gt.nres
- & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
- call ssbond_ene(iii,jjj,eij)
- ehpb=ehpb+2*eij
+ 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+eij
endif
cd write (iout,*) "eij",eij
- else
+cd & ' waga=',waga,' fac=',fac
+! else if (ii.gt.nres .and. jj.gt.nres) then
+ else
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
+ if (irestr_type(i).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,6f10.3,i5)')
+ & "edisL",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+ & ehpb,irestr_type(i)
+ else if (irestr_type(i).eq.10) then
+c AL 6//19/2018 cross-link restraints
+ xdis = 0.5d0*(dd/forcon(i))**2
+ expdis = dexp(-xdis)
+c aux=(dhpb(i)+dhpb1(i)*xdis)*expdis+fordepth(i)
+ aux=(dhpb(i)+dhpb1(i)*xdis*xdis)*expdis+fordepth(i)
+c write (iout,*)"HERE: xdis",xdis," expdis",expdis," aux",aux,
+c & " wboltzd",wboltzd
+ ehpb=ehpb-wboltzd*xlscore(i)*dlog(aux)
+c fac=-wboltzd*(dhpb1(i)*(1.0d0-xdis)-dhpb(i))
+ fac=-wboltzd*xlscore(i)*(dhpb1(i)*(2.0d0-xdis)*xdis-dhpb(i))
+ & *expdis/(aux*forcon(i)**2)
+ if (energy_dec) write(iout,'(a6,2i5,6f10.3,i5)')
+ & "edisX",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+ & -wboltzd*xlscore(i)*dlog(aux),irestr_type(i)
+ else if (irestr_type(i).eq.2) then
+c Quartic restraints
+ ehpb=ehpb+forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)')
+ & "edisQ",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+ & forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)),irestr_type(i)
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+ else
+c Quadratic restraints
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
+ ehpb=ehpb+0.5d0*waga*rdis*rdis
+ if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)')
+ & "edisS",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+ & 0.5d0*waga*rdis*rdis,irestr_type(i)
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
+ endif
+c Calculate Cartesian gradient
+ 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
- do j=1,3
- ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
- ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
- enddo
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
endif
cgrad do j=iii,jjj-1
cgrad do k=1,3
enddo
endif
enddo
- ehpb=0.5D0*ehpb
return
end
C--------------------------------------------------------------------------
c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
c endif
enddo
+
estr=0.5d0*AKP*estr+estr1
c
c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
enddo
+
C Ufff.... We've done all this!!!
return
end
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)
enddo
+
return
end
#endif
c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' energy=',ekont*ees,
c & 'gradcorr_long'
C Calculate the multi-body contribution to energy.
-c ecorr=ecorr+ekont*ees
+C ecorr=ecorr+ekont*ees
C Calculate multi-body contributions to the gradient.
coeffpees0pij=coeffp*ees0pij
coeffmees0mij=coeffm*ees0mij
else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then
gamm=(rr-(r_cut-rlamb))/rlamb
sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0)
- else
+ else
sscale2=0d0
endif
- return
- end
+ return
+ end
C-----------------------------------------------------------------------
double precision function sscalgrad2(r,r_cut,r0,rlamb)
implicit none
else
sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb
endif
- else
+ else
sscalgrad2=0.0d0
endif
- return
- end
+ return
+ end
c----------------------------------------------------------------------------
subroutine e_saxs(Esaxs_constr)
implicit none
enddo ! i
#ifdef MPI
if (nfgtasks.gt.1) then
- call MPI_Reduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION,
- & MPI_SUM,king,FG_COMM,IERR)
- if (fg_rank.eq.king) then
+ call MPI_AllReduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION,
+ & MPI_SUM,FG_COMM,IERR)
+c if (fg_rank.eq.king) then
do k=1,nsaxs
Pcalc(k) = Pcalc_(k)
enddo
- endif
- call MPI_Reduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres,
- & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
- if (fg_rank.eq.king) then
- do i=1,nres
- do l=1,3
- do k=1,nsaxs
- PgradC(k,l,i) = PgradC_(k,l,i)
- enddo
- enddo
- enddo
- endif
+c endif
+c call MPI_AllReduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres,
+c & MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
+c if (fg_rank.eq.king) then
+c do i=1,nres
+c do l=1,3
+c do k=1,nsaxs
+c PgradC(k,l,i) = PgradC_(k,l,i)
+c enddo
+c enddo
+c enddo
+c endif
#ifdef ALLSAXS
- call MPI_Reduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres,
- & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
- if (fg_rank.eq.king) then
- do i=1,nres
- do l=1,3
- do k=1,nsaxs
- PgradX(k,l,i) = PgradX_(k,l,i)
- enddo
- enddo
- enddo
- endif
+c call MPI_AllReduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres,
+c & MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
+c if (fg_rank.eq.king) then
+c do i=1,nres
+c do l=1,3
+c do k=1,nsaxs
+c PgradX(k,l,i) = PgradX_(k,l,i)
+c enddo
+c enddo
+c enddo
+c endif
#endif
endif
#endif
-#ifdef MPI
- if (fg_rank.eq.king) then
-#endif
Cnorm = 0.0d0
do k=1,nsaxs
Cnorm = Cnorm + Pcalc(k)
enddo
+#ifdef MPI
+ if (fg_rank.eq.king) then
+#endif
Esaxs_constr = dlog(Cnorm)-wsaxs0
do k=1,nsaxs
if (Pcalc(k).gt.0.0d0)
#ifdef DEBUG
write (iout,*) "Cnorm",Cnorm," Esaxs_constr",Esaxs_constr
#endif
+#ifdef MPI
+ endif
+#endif
+ gsaxsC=0.0d0
+ gsaxsX=0.0d0
do i=nnt,nct
do l=1,3
auxC=0.0d0
enddo
enddo
#ifdef MPI
- endif
+c endif
#endif
return
end
projects / unres.git / blobdiff