weights_(17)=wbond
weights_(18)=scal14
weights_(21)=wsccor
+ weights_(25)=wsaxs
C FG Master broadcasts the WEIGHTS_ array
call MPI_Bcast(weights_(1),n_ene,
& MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
wbond=weights(17)
scal14=weights(18)
wsccor=weights(21)
+ wsaxs=weights(25)
endif
time_Bcast=time_Bcast+MPI_Wtime()-time00
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
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
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
cd write (iout,*) "multibody_hb ecorr",ecorr
endif
+c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode
+ if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
+ call e_saxs(Esaxs_constr)
+c write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr
+ else if (nsaxs.gt.0 .and. saxs_mode.gt.0) then
+ call e_saxsC(Esaxs_constr)
+c write (iout,*) "From EsaxsC: Esaxs_constr",Esaxs_constr
+ else
+ Esaxs_constr = 0.0d0
+ 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
energia(22)=eliptran
energia(23)=Eafmforce
energia(24)=ehomology_constr
+ energia(25)=Esaxs_constr
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"
eliptran=energia(22)
Eafmforce=energia(23)
ehomology_constr=energia(24)
+ esaxs_constr=energia(25)
+c write (iout,*) "sum_energy esaxs_constr",esaxs_constr,
+c & " wsaxs",wsaxs
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
- & +wliptran*eliptran+Eafmforce
+ & +wsaxs*esaxs_constr+wliptran*eliptran+Eafmforce
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
- & +wliptran*eliptran
+ & +wsaxs*esaxs_constr+wliptran*eliptran
& +Eafmforce
#endif
energia(0)=etot
#endif
#ifdef DEBUG
write (iout,*) "sum_gradient gvdwc, gvdwx"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)')
+ do i=0,nres
+ write (iout,'(i3,3e15.5,5x,3e15.5)')
& i,(gvdwx(j,i),j=1,3),(gvdwc(j,i),j=1,3)
enddo
call flush(iout)
#endif
+#ifdef DEBUG
+ write (iout,*) "sum_gradient gsaxsc, gsaxsx"
+ do i=0,nres
+ 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)
+#endif
#ifdef MPI
C FG slaves call the following matching MPI_Bcast in ERGASTULUM
if (nfgtasks.gt.1 .and. fg_rank.eq.0)
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wsaxs*gsaxsc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wsaxs*gsaxsc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,-1,-1
+ do i=nres-2,0,-1
+c do i=nres-2,nnt,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
#ifdef DEBUG
write (iout,*) "gradbufc after summing"
do i=1,nres
- write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+ write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3)
enddo
call flush(iout)
#endif
#endif
#ifdef DEBUG
write (iout,*) "gradbufc"
- do i=1,nres
- write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+ do i=0,nres
+ write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3)
enddo
call flush(iout)
#endif
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,-1,-1
+ do i=nres-2,0,-1
+c do i=nres-2,nnt,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
c enddo
#ifdef DEBUG
write (iout,*) "gradbufc after summing"
- do i=1,nres
- write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+ do i=0,nres
+ write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3)
enddo
call flush(iout)
#endif
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)+
#endif
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
- & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
- & wsccor*gsccorx(j,i)
+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
+ & +wsaxs*gsaxsx(j,i)
+ & +wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
& +wliptran*gliptranx(j,i)
enddo
#ifdef DEBUG
write (iout,*) "gradc gradx gloc"
do i=1,nres
- write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)')
+ write (iout,'(i5,3e15.5,5x,3e15.5,5x,e15.5)')
& i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg)
enddo
#endif
Uconst=energia(20)
esccor=energia(21)
ehomology_constr=energia(24)
+ esaxs_constr=energia(25)
eliptran=energia(22)
Eafmforce=energia(23)
#ifdef SPLITELE
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
- & edihcnstr,ehomology_constr, ebr*nss,
+ & edihcnstr,ehomology_constr,esaxs_constr*wsaxs, ebr*nss,
& Uconst,eliptran,wliptran,Eafmforce,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
+ & 'E_SAXS=',1pE16.6,' (SAXS restraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ehomology_constr,ebr*nss,Uconst,
+ & ehomology_constr,esaxs_constr*wsaxs,ebr*nss,Uconst,
& eliptran,wliptran,Eafmforc,
& etot
10 format (/'Virtual-chain energies:'//
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
+ & 'E_SAXS=',1pE16.6,' (SAXS restraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
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,*) '*******i',i,' iti1',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)
c if (i .gt. iatel_s+2) then
if (i .gt. nnt+2) then
mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
-c write (iout,*) 'mu ',mu(:,i-2),i-2
+cd write (iout,*) 'mu ',mu(:,i-2),i-2
+cd write (iout,*) 'b1 ',b1(:,i-1),i-2
+cd write (iout,*) 'Ub2 ',Ub2(:,i-2),i-2
+cd write (iout,*) 'Ug ',Ug(:,:,i-2),i-2
+cd write (iout,*) 'b2 ',b2(:,i-2),i-2
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
C
C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
do i=iturn3_start,iturn3_end
- if (i.le.1) cycle
+CAna if (i.le.1) cycle
C write(iout,*) "tu jest i",i
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
- & .or.((i+4).gt.nres)
- & .or.((i-1).le.0)
+CAna & .or.((i+4).gt.nres)
+CAna & .or.((i-1).le.0)
C end of changes by Ana
& .or. itype(i+2).eq.ntyp1
& .or. itype(i+3).eq.ntyp1) cycle
- if(i.gt.1)then
- if(itype(i-1).eq.ntyp1)cycle
- end if
- if(i.LT.nres-3)then
- if (itype(i+4).eq.ntyp1) cycle
- end if
+CAna if(i.gt.1)then
+CAna if(itype(i-1).eq.ntyp1)cycle
+CAna end if
+CAna if(i.LT.nres-3)then
+CAna if (itype(i+4).eq.ntyp1) cycle
+CAna end if
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
- if (i.le.1) cycle
+cAna if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
- & .or.((i+5).gt.nres)
- & .or.((i-1).le.0)
+cAna & .or.((i+5).gt.nres)
+cAna & .or.((i-1).le.0)
C end of changes suggested by Ana
& .or. itype(i+3).eq.ntyp1
& .or. itype(i+4).eq.ntyp1
- & .or. itype(i+5).eq.ntyp1
- & .or. itype(i).eq.ntyp1
- & .or. itype(i-1).eq.ntyp1
+cAna & .or. itype(i+5).eq.ntyp1
+cAna & .or. itype(i).eq.ntyp1
+cAna & .or. itype(i-1).eq.ntyp1
& ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
- if (i.le.1) cycle
+cAna if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
- & .or.((i+2).gt.nres)
- & .or.((i-1).le.0)
+cAna & .or.((i+2).gt.nres)
+cAna & .or.((i-1).le.0)
C end of changes by Ana
- & .or. itype(i+2).eq.ntyp1
- & .or. itype(i-1).eq.ntyp1
+cAna & .or. itype(i+2).eq.ntyp1
+cAna & .or. itype(i-1).eq.ntyp1
& ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
C write (iout,*) i,j
- if (j.le.1) cycle
+cAna if (j.le.1) cycle
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
- & .or.((j+2).gt.nres)
- & .or.((j-1).le.0)
+cAna & .or.((j+2).gt.nres)
+cAna & .or.((j-1).le.0)
C end of changes by Ana
- & .or.itype(j+2).eq.ntyp1
- & .or.itype(j-1).eq.ntyp1
+cAna & .or.itype(j+2).eq.ntyp1
+cAna & .or.itype(j-1).eq.ntyp1
&) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
+
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
C Cartesian derivatives
+!DIR$ UNROLL(0)
do l=1,3
c ghalf1=0.5d0*agg(l,1)
c ghalf2=0.5d0*agg(l,2)
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
+ 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
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+2*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
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
- if (i.eq.2) cycle
+c if (i.eq.2) cycle
c print *,i,itype(i-1),itype(i),itype(i-2)
if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1)
& .or.(itype(i).eq.ntyp1)) cycle
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
enddo
else
phii=0.0d0
- ityp1=nthetyp+1
+ ityp1=ithetyp(itype(i-2))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp(itype(i))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
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 & dscp1,dscp2,sumene
c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+ & 'escloc',i,sumene
c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
c & ,zz,xx,yy
c#define DEBUG
c
- do i=1,19
+ do i=1,max_template
distancek(i)=9999999.9
enddo
j = jres_homo(ii)
dij=dist(i,j)
c write (iout,*) "dij(",i,j,") =",dij
+ nexl=0
do k=1,constr_homology
c write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii)
- if(.not.l_homo(k,ii)) cycle
+ if(.not.l_homo(k,ii)) then
+ nexl=nexl+1
+ cycle
+ endif
distance(k)=odl(k,ii)-dij
c write (iout,*) "distance(",k,") =",distance(k)
c
endif
enddo
- min_odl=minval(distancek)
+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)
+ & min_odl=distancek(kk)
+ enddo
+
c write (iout,* )"min_odl",min_odl
#ifdef DEBUG
write (iout,*) "ij dij",i,j,dij
write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
write (iout,* )"min_odl",min_odl
#endif
+#ifdef OLDRESTR
odleg2=0.0d0
+#else
+ if (waga_dist.ge.0.0d0) then
+ odleg2=nexl
+ else
+ odleg2=0.0d0
+ endif
+#endif
do k=1,constr_homology
c Nie wiem po co to liczycie jeszcze raz!
c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
do i=idihconstr_start_homo,idihconstr_end_homo
kat2=0.0d0
c betai=beta(i,i+1,i+2,i+3)
- betai = phi(i+3)
+ betai = phi(i)
c write (iout,*) "betai =",betai
do k=1,constr_homology
dih_diff(k)=pinorm(dih(k,i)-betai)
-c write (iout,*) "dih_diff(",k,") =",dih_diff(k)
+cd write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k)
+cd & ,sigma_dih(k,i)
c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
c & -(6.28318-dih_diff(i,k))
c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
c & 6.28318+dih_diff(i,k)
-
+#ifdef OLD_DIHED
kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#else
+ kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#endif
c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
gdih(k)=dexp(kat3)
kat2=kat2+gdih(k)
sum_gdih=kat2
sum_sgdih=0.0d0
do k=1,constr_homology
+#ifdef OLD_DIHED
sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+#else
+ sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd
+#endif
c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
sum_sgdih=sum_sgdih+sgdih
enddo
c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
ccc & gloc(nphi+i-3,icg)
- gloc(i,icg)=gloc(i,icg)+grad_dih3
+ gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3
c if (i.eq.25) then
c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
c endif
c dtheta_i=theta(j)-thetaref(j,iref)
c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing
theta_diff(k)=thetatpl(k,i)-theta(i)
+cd write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k)
+cd & ,sigma_theta(k,i)
+
c
utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
- gutheta_i=gutheta_i+dexp(utheta_i) ! Sum of Gaussians (pk)
+ gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk)
c Gradient for single Gaussian restraint in subr Econstr_back
c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
c
dyy=-yytpl(k,i)+yytab(i) ! ibid y
dzz=-zztpl(k,i)+zztab(i) ! ibid z
c write(iout,*) "dxx, dyy, dzz"
-c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+cd write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i)
c
usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument
c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
c uscdiffk(k)=usc_diff(i)
guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
- guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk)
+c write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i),
+c & " guscdiff2",guscdiff2(k)
+ guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk)
c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
c & xxref(j),yyref(j),zzref(j)
enddo
esccor=0.0D0
do i=itau_start,itau_end
if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
- esccor_ii=0.0D0
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
c write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
phii=phi(i)
do intertyp=1,3 !intertyp
+ esccor_ii=0.0D0
cc Added 09 May 2012 (Adasko)
cc Intertyp means interaction type of backbone mainchain correlation:
c 1 = SC...Ca...Ca...Ca
v2ij=v2sccor(j,intertyp,isccori,isccori1)
cosphi=dcos(j*tauangle(intertyp,i))
sinphi=dsin(j*tauangle(intertyp,i))
+ if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
+ if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)')
+ & 'esccor',i,intertyp,esccor_ii
c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
if (lprn)
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
return
end
+c----------------------------------------------------------------------------
+ double precision function sscale2(r,r_cut,r0,rlamb)
+ implicit none
+ double precision r,gamm,r_cut,r0,rlamb,rr
+ rr = dabs(r-r0)
+c write (2,*) "r",r," r_cut",r_cut," r0",r0," rlamb",rlamb
+c write (2,*) "rr",rr
+ if(rr.lt.r_cut-rlamb) then
+ sscale2=1.0d0
+ 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
+ sscale2=0d0
+ endif
+ return
+ end
+C-----------------------------------------------------------------------
+ double precision function sscalgrad2(r,r_cut,r0,rlamb)
+ implicit none
+ double precision r,gamm,r_cut,r0,rlamb,rr
+ rr = dabs(r-r0)
+ if(rr.lt.r_cut-rlamb) then
+ sscalgrad2=0.0d0
+ else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then
+ gamm=(rr-(r_cut-rlamb))/rlamb
+ if (r.ge.r0) then
+ sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb
+ endif
+ else
+ sscalgrad2=0.0d0
+ endif
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine e_saxs(Esaxs_constr)
+ implicit none
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+ include "COMMON.SETUP"
+ integer IERR
+#endif
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.GEO'
+ include 'COMMON.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.CONTROL'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ include 'COMMON.FFIELD'
+c
+ double precision Esaxs_constr
+ integer i,iint,j,k,l
+ double precision PgradC(maxSAXS,3,maxres),
+ & PgradX(maxSAXS,3,maxres),Pcalc(maxSAXS)
+#ifdef MPI
+ double precision PgradC_(maxSAXS,3,maxres),
+ & PgradX_(maxSAXS,3,maxres),Pcalc_(maxSAXS)
+#endif
+ double precision dk,dijCACA,dijCASC,dijSCCA,dijSCSC,
+ & sigma2CACA,sigma2CASC,sigma2SCCA,sigma2SCSC,expCACA,expCASC,
+ & expSCCA,expSCSC,CASCgrad,SCCAgrad,SCSCgrad,aux,auxC,auxC1,
+ & auxX,auxX1,CACAgrad,Cnorm
+ double precision sss2,ssgrad2,rrr,sscalgrad2,sscale2
+ double precision dist
+ external dist
+c SAXS restraint penalty function
+#ifdef DEBUG
+ write(iout,*) "------- SAXS penalty function start -------"
+ write (iout,*) "nsaxs",nsaxs
+ write (iout,*) "Esaxs: iatsc_s",iatsc_s," iatsc_e",iatsc_e
+ write (iout,*) "Psaxs"
+ do i=1,nsaxs
+ write (iout,'(i5,e15.5)') i, Psaxs(i)
+ enddo
+#endif
+ Esaxs_constr = 0.0d0
+ do k=1,nsaxs
+ Pcalc(k)=0.0d0
+ do j=1,nres
+ do l=1,3
+ PgradC(k,l,j)=0.0d0
+ PgradX(k,l,j)=0.0d0
+ enddo
+ enddo
+ enddo
+ do i=iatsc_s,iatsc_e
+ if (itype(i).eq.ntyp1) cycle
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ if (itype(j).eq.ntyp1) cycle
+#ifdef ALLSAXS
+ dijCACA=dist(i,j)
+ dijCASC=dist(i,j+nres)
+ dijSCCA=dist(i+nres,j)
+ dijSCSC=dist(i+nres,j+nres)
+ sigma2CACA=2.0d0/(pstok**2)
+ sigma2CASC=4.0d0/(pstok**2+restok(itype(j))**2)
+ sigma2SCCA=4.0d0/(pstok**2+restok(itype(i))**2)
+ sigma2SCSC=4.0d0/(restok(itype(j))**2+restok(itype(i))**2)
+ do k=1,nsaxs
+ dk = distsaxs(k)
+ expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+ if (itype(j).ne.10) then
+ expCASC = dexp(-0.5d0*sigma2CASC*(dijCASC-dk)**2)
+ else
+ endif
+ expCASC = 0.0d0
+ if (itype(i).ne.10) then
+ expSCCA = dexp(-0.5d0*sigma2SCCA*(dijSCCA-dk)**2)
+ else
+ expSCCA = 0.0d0
+ endif
+ if (itype(i).ne.10 .and. itype(j).ne.10) then
+ expSCSC = dexp(-0.5d0*sigma2SCSC*(dijSCSC-dk)**2)
+ else
+ expSCSC = 0.0d0
+ endif
+ Pcalc(k) = Pcalc(k)+expCACA+expCASC+expSCCA+expSCSC
+#ifdef DEBUG
+ write(iout,*) "i j k Pcalc",i,j,Pcalc(k)
+#endif
+ CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+ CASCgrad = sigma2CASC*(dijCASC-dk)*expCASC
+ SCCAgrad = sigma2SCCA*(dijSCCA-dk)*expSCCA
+ SCSCgrad = sigma2SCSC*(dijSCSC-dk)*expSCSC
+ do l=1,3
+c CA CA
+ aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+c CA SC
+ if (itype(j).ne.10) then
+ aux = CASCgrad*(C(l,j+nres)-C(l,i))/dijCASC
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+ PgradX(k,l,j) = PgradX(k,l,j)+aux
+ endif
+c SC CA
+ if (itype(i).ne.10) then
+ aux = SCCAgrad*(C(l,j)-C(l,i+nres))/dijSCCA
+ PgradX(k,l,i) = PgradX(k,l,i)-aux
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+ endif
+c SC SC
+ if (itype(i).ne.10 .and. itype(j).ne.10) then
+ aux = SCSCgrad*(C(l,j+nres)-C(l,i+nres))/dijSCSC
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+ PgradX(k,l,i) = PgradX(k,l,i)-aux
+ PgradX(k,l,j) = PgradX(k,l,j)+aux
+ endif
+ enddo ! l
+ enddo ! k
+#else
+ dijCACA=dist(i,j)
+ sigma2CACA=scal_rad**2*0.25d0/
+ & (restok(itype(j))**2+restok(itype(i))**2)
+
+ IF (saxs_cutoff.eq.0) THEN
+ do k=1,nsaxs
+ dk = distsaxs(k)
+ expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+ Pcalc(k) = Pcalc(k)+expCACA
+ CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+ do l=1,3
+ aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+ enddo ! l
+ enddo ! k
+ ELSE
+ rrr = saxs_cutoff*2.0d0/dsqrt(sigma2CACA)
+ do k=1,nsaxs
+ dk = distsaxs(k)
+c write (2,*) "ijk",i,j,k
+ sss2 = sscale2(dijCACA,rrr,dk,0.3d0)
+ if (sss2.eq.0.0d0) cycle
+ ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0)
+ if (energy_dec) write(iout,'(a4,3i5,5f10.4)')
+ & 'saxs',i,j,k,dijCACA,rrr,dk,sss2,ssgrad2
+ expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)*sss2
+ Pcalc(k) = Pcalc(k)+expCACA
+#ifdef DEBUG
+ write(iout,*) "i j k Pcalc",i,j,Pcalc(k)
+#endif
+ CACAgrad = -sigma2CACA*(dijCACA-dk)*expCACA+
+ & ssgrad2*expCACA/sss2
+ do l=1,3
+c CA CA
+ aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA
+ PgradC(k,l,i) = PgradC(k,l,i)+aux
+ PgradC(k,l,j) = PgradC(k,l,j)-aux
+ enddo ! l
+ enddo ! k
+ ENDIF
+#endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+#ifdef MPI
+ if (nfgtasks.gt.1) 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
+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
+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
+ 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)
+ & Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k))
+#ifdef DEBUG
+ write (iout,*) "k",k," Esaxs_constr",Esaxs_constr
+#endif
+ enddo
+#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
+ auxC1=0.0d0
+ auxX=0.0d0
+ auxX1=0.d0
+ do k=1,nsaxs
+ if (Pcalc(k).gt.0)
+ & auxC = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k)
+ auxC1 = auxC1+PgradC(k,l,i)
+#ifdef ALLSAXS
+ auxX = auxX +Psaxs(k)*PgradX(k,l,i)/Pcalc(k)
+ auxX1 = auxX1+PgradX(k,l,i)
+#endif
+ enddo
+ gsaxsC(l,i) = auxC - auxC1/Cnorm
+#ifdef ALLSAXS
+ gsaxsX(l,i) = auxX - auxX1/Cnorm
+#endif
+c write (iout,*) "l i",l,i," gradC",wsaxs*(auxC - auxC1/Cnorm),
+c * " gradX",wsaxs*(auxX - auxX1/Cnorm)
+ enddo
+ enddo
+#ifdef MPI
+c endif
+#endif
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine e_saxsC(Esaxs_constr)
+ implicit none
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+ include "COMMON.SETUP"
+ integer IERR
+#endif
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.GEO'
+ include 'COMMON.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.CONTROL'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ include 'COMMON.FFIELD'
+c
+ double precision Esaxs_constr
+ integer i,iint,j,k,l
+ double precision PgradC(3,maxres),PgradX(3,maxres),Pcalc,logPtot
+#ifdef MPI
+ double precision gsaxsc_(3,maxres),gsaxsx_(3,maxres),logPtot_
+#endif
+ double precision dk,dijCASPH,dijSCSPH,
+ & sigma2CA,sigma2SC,expCASPH,expSCSPH,
+ & CASPHgrad,SCSPHgrad,aux,auxC,auxC1,
+ & auxX,auxX1,Cnorm
+c SAXS restraint penalty function
+#ifdef DEBUG
+ write(iout,*) "------- SAXS penalty function start -------"
+ write (iout,*) "nsaxs",nsaxs
+
+ do i=nnt,nct
+ print *,MyRank,"C",i,(C(j,i),j=1,3)
+ enddo
+ do i=nnt,nct
+ print *,MyRank,"CSaxs",i,(Csaxs(j,i),j=1,3)
+ enddo
+#endif
+ Esaxs_constr = 0.0d0
+ logPtot=0.0d0
+ do j=isaxs_start,isaxs_end
+ Pcalc=0.0d0
+ do i=1,nres
+ do l=1,3
+ PgradC(l,i)=0.0d0
+ PgradX(l,i)=0.0d0
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).eq.ntyp1) cycle
+ dijCASPH=0.0d0
+ dijSCSPH=0.0d0
+ do l=1,3
+ dijCASPH=dijCASPH+(C(l,i)-Csaxs(l,j))**2
+ enddo
+ if (itype(i).ne.10) then
+ do l=1,3
+ dijSCSPH=dijSCSPH+(C(l,i+nres)-Csaxs(l,j))**2
+ enddo
+ endif
+ sigma2CA=2.0d0/pstok**2
+ sigma2SC=4.0d0/restok(itype(i))**2
+ expCASPH = dexp(-0.5d0*sigma2CA*dijCASPH)
+ expSCSPH = dexp(-0.5d0*sigma2SC*dijSCSPH)
+ Pcalc = Pcalc+expCASPH+expSCSPH
+#ifdef DEBUG
+ write(*,*) "processor i j Pcalc",
+ & MyRank,i,j,dijCASPH,dijSCSPH, Pcalc
+#endif
+ CASPHgrad = sigma2CA*expCASPH
+ SCSPHgrad = sigma2SC*expSCSPH
+ do l=1,3
+ aux = (C(l,i+nres)-Csaxs(l,j))*SCSPHgrad
+ PgradX(l,i) = PgradX(l,i) + aux
+ PgradC(l,i) = PgradC(l,i)+(C(l,i)-Csaxs(l,j))*CASPHgrad+aux
+ enddo ! l
+ enddo ! i
+ do i=nnt,nct
+ do l=1,3
+ gsaxsc(l,i)=gsaxsc(l,i)+PgradC(l,i)/Pcalc
+ gsaxsx(l,i)=gsaxsx(l,i)+PgradX(l,i)/Pcalc
+ enddo
+ enddo
+ logPtot = logPtot - dlog(Pcalc)
+c print *,"me",me,MyRank," j",j," logPcalc",-dlog(Pcalc),
+c & " logPtot",logPtot
+ enddo ! j
+#ifdef MPI
+ if (nfgtasks.gt.1) then
+c write (iout,*) "logPtot before reduction",logPtot
+ call MPI_Reduce(logPtot,logPtot_,1,MPI_DOUBLE_PRECISION,
+ & MPI_SUM,king,FG_COMM,IERR)
+ logPtot = logPtot_
+c write (iout,*) "logPtot after reduction",logPtot
+ call MPI_Reduce(gsaxsC(1,1),gsaxsC_(1,1),3*nres,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ if (fg_rank.eq.king) then
+ do i=1,nres
+ do l=1,3
+ gsaxsC(l,i) = gsaxsC_(l,i)
+ enddo
+ enddo
+ endif
+ call MPI_Reduce(gsaxsX(1,1),gsaxsX_(1,1),3*nres,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ if (fg_rank.eq.king) then
+ do i=1,nres
+ do l=1,3
+ gsaxsX(l,i) = gsaxsX_(l,i)
+ enddo
+ enddo
+ endif
+ endif
+#endif
+ Esaxs_constr = logPtot
+ return
+ end
projects / unres.git / blobdiff