enddo
endif
+ if(homol_nset.gt.1) then
+ i_econstr=24
+ nset=homol_nset
+ do i=1,nset
+ mset(i)=1
+ enddo
+ endif
+
+ if(usampl) i_econstr=20
+
k=0
rep2i(k)=-1
do il=1,max0(nset,1)
read (irest2,*) ndowna(0,il),
& (ndowna(i,il),i=1,ndowna(0,il))
enddo
- if(usampl.or.hremd.gt.0) then
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
read (irest2,*)
read (irest2,*) nset
read (irest2,*)
if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
endif
- if(usampl.or.hremd.gt.0) then
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
iset=i2set(me)
if (hremd.gt.0) call set_hweights(iset)
if(me.eq.king.or..not.out1file)
write (irest1,*) ndowna(0,il),
& (ndowna(i,il),i=1,ndowna(0,il))
enddo
- if(usampl.or.hremd.gt.0) then
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
write (irest1,*) "nset"
write (irest1,*) nset
write (irest1,*) "mset"
do i=1,2*nres
write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
- if(usampl.or.hremd.gt.0) then
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
write (irest2,*) iset
endif
close(irest2)
c & CG_COMM,ierr)
potEcomp(n_ene+1)=t_bath
t_bath_old=t_bath
- if (usampl) then
+ if (usampl.or.homol_nset.gt.1) then
potEcomp(n_ene+2)=iset
+c write(iout,*) potEcomp(n_ene+1),potEcomp(n_ene+2),iset,nset
if (iset.lt.nset) then
i_set_temp=iset
iset=iset+1
- call EconstrQ
- potEcomp(n_ene+3)=Uconst
+ if (homol_nset.gt.1) then
+ call e_modeller(potEcomp(n_ene+3))
+c write(iout,*) "iset+1",potEcomp(n_ene+3)
+ else
+ call EconstrQ
+ potEcomp(n_ene+3)=Uconst
+ endif
iset=i_set_temp
endif
if (iset.gt.1) then
i_set_temp=iset
iset=iset-1
- call EconstrQ
- potEcomp(n_ene+4)=Uconst
+ if (homol_nset.gt.1) then
+ call e_modeller(potEcomp(n_ene+4))
+c write(iout,*) "iset-1",potEcomp(n_ene+4)
+ else
+ call EconstrQ
+ potEcomp(n_ene+4)=Uconst
+ endif
iset=i_set_temp
endif
endif
if (me.eq.king) then
+ if(homol_nset.gt.1) write(iout,*)
+ & 'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
do i=1,nodes
remd_t_bath(i)=remd_ene(n_ene+1,i)
iremd_iset(i)=remd_ene(n_ene+2,i)
+ if(homol_nset.gt.1)
+ & write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
+ & i,remd_ene(i_econstr,i),
+ & remd_ene(n_ene+1,i),iremd_iset(i),
+ & remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
enddo
#ifdef DEBUG
if(lmuca) then
endif
#endif
c-------------------------------------
- IF(.not.usampl.and.hremd.eq.0) THEN
+ IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
#ifdef DEBUG
write (iout,*) "Enter exchnge, remd_m",remd_m(1),
& " nodes",nodes
enddo
cd write (iout,*) "exchange completed"
cd call flush(iout)
- ELSEIF (usampl) THEN
+ ELSEIF (usampl.or.homol_nset.gt.1) THEN
do ii=1,nodes
-cd write(iout,*) "########",ii
+c write(iout,*) "########",ii
i_temp=iran_num(1,nrep)
i_mult=iran_num(1,remd_m(i_temp))
i_mset=iran_num(1,mset(i_iset))
i=i_index(i_temp,i_mult,i_iset,i_mset)
-cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
+c write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
if(t_exchange_only)then
i_dir=1
else
i_dir=iran_num(1,3)
endif
-cd write(iout,*) "i_dir=",i_dir
+c write(iout,*) "i_dir=",i_dir
if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
i_iset1=i_iset+1
i_mset1=iran_num(1,mset(i_iset1))
iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
- econstr_temp_i=remd_ene(20,i)
- econstr_temp_iex=remd_ene(20,iex)
- remd_ene(20,i)=remd_ene(n_ene+3,i)
- remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+
+ econstr_temp_i=remd_ene(i_econstr,i)
+ econstr_temp_iex=remd_ene(i_econstr,iex)
+ remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+ remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
i_iset1=i_iset+1
i_mset1=iran_num(1,mset(i_iset1))
iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
- econstr_temp_i=remd_ene(20,i)
- econstr_temp_iex=remd_ene(20,iex)
- remd_ene(20,i)=remd_ene(n_ene+3,i)
- remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+ econstr_temp_i=remd_ene(i_econstr,i)
+ econstr_temp_iex=remd_ene(i_econstr,iex)
+ remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+ remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
else
goto 444
endif
-cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
+c write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
ctime call flush(iout)
c Swap temperatures between conformations i and iex with recalculating the free energies
call sum_energy(remd_ene(0,iex),.false.)
ene_iex_i=remd_ene(0,iex)
-cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
+cdebug
+c ERROR only makes sense for dir =1
+c write (iout,*) "ene_iex_i",remd_ene(0,iex)
c call sum_energy(remd_ene(0,i),.false.)
-cd write (iout,*) "ene_i_i",remd_ene(0,i)
+c write (iout,*) "ene_i_i",remd_ene(0,i)
c write (iout,*) "rescaling weights with temperature",
c & remd_t_bath(iex)
c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
-c write (iout,*) "ERROR: inconsistent energies:",i,
+c write (iout,*) "ERROR: inconsistent energies i:",i,
c & ene_i_i,remd_ene(0,i)
c endif
+cdebug_end
call rescale_weights(remd_t_bath(iex))
call sum_energy(remd_ene(0,i),.false.)
cd write (iout,*) "ene_i_iex",remd_ene(0,i)
ene_i_iex=remd_ene(0,i)
+cdebug
+c ERROR only makes sense for dir =1
c call sum_energy(remd_ene(0,iex),.false.)
c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
-c write (iout,*) "ERROR: inconsistent energies:",iex,
+c write (iout,*) "ERROR: inconsistent energies iex:",iex,
c & ene_iex_iex,remd_ene(0,iex)
c endif
-cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
+c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
c write (iout,*) "i",i," iex",iex
-cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-cd & " ene_i_iex",ene_i_iex,
-cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
+c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
+c & " ene_i_iex",ene_i_iex,
+c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
+cdebug_end
delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
& (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
delta=-delta
-cd write(iout,*) 'delta',delta
+c write(iout,*) 'delta',delta
c delta=(remd_t_bath(i)-remd_t_bath(iex))*
c & (remd_ene(i)-remd_ene(iex))/Rb/
c & (remd_t_bath(i)*remd_t_bath(iex))
& iremd_tot_usa(int(i2set(i-1)))=
& iremd_tot_usa(int(i2set(i-1)))+1
xxx=ran_number(0.0d0,1.0d0)
-cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
+c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
if (delta .gt. xxx) then
tmp=remd_t_bath(i)
remd_t_bath(i)=remd_t_bath(iex)
else
remd_ene(0,iex)=ene_iex_iex
remd_ene(0,i)=ene_i_i
- remd_ene(20,iex)=econstr_temp_iex
- remd_ene(20,i)=econstr_temp_i
+ remd_ene(i_econstr,iex)=econstr_temp_iex
+ remd_ene(i_econstr,i)=econstr_temp_i
endif
cd do il=1,nset
& ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
do i=1,nset
if(iremd_tot_usa(i).ne.0)
& write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
& CG_COMM,ierr)
cd write (iout,*) "After scatter"
cd call flush(iout)
- if(usampl.or.hremd.gt.0)
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1)
& call mpi_scatter(iremd_iset,1,mpi_integer,
& iset,1,mpi_integer,king,
& CG_COMM,ierr)
include 'COMMON.CHAIN'
include 'COMMON.SBRIDGE'
include 'COMMON.INTERACT'
+ include 'COMMON.CONTROL'
real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
& d_restart2(3,2*maxres*maxprocs)
enddo
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
call xdrfint_(ixdrf, nset, iret)
do i=1,nset
call xdrfint_(ixdrf,mset(i), iret)
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
call xdrfint(ixdrf, nset, iret)
do i=1,nset
call xdrfint(ixdrf,mset(i), iret)
call xdrfint_(ixdrf, nss, iret)
do j=1,nss
if (dyn_ss) then
- call xdrfint(ixdrf, idssb(j)+nres, iret)
- call xdrfint(ixdrf, jdssb(j)+nres, iret)
+ call xdrfint_(ixdrf, idssb(j)+nres, iret)
+ call xdrfint_(ixdrf, jdssb(j)+nres, iret)
else
call xdrfint_(ixdrf, ihpb(j), iret)
call xdrfint_(ixdrf, jhpb(j), iret)
include 'COMMON.CHAIN'
include 'COMMON.SBRIDGE'
include 'COMMON.INTERACT'
+ include 'COMMON.CONTROL'
real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
& t5_restart1(5)
integer*2 i_index
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
#ifdef AIX
if(me.eq.king)then
call xdrfint_(ixdrf, nset, iret)
enddo
endif
#endif
++c<<<<<<< HEAD
+c Corrected AL 8/19/2014: each processor needs whole iset array not only its
++cv=======
+ Corrected AL 8/19/2014: each processor needs whole iset array not only its
++c>>>>>>> 66693b0684c228404e7aadcffe6d2a8c9f489063
c own element
c call mpi_scatter(i2set,1,mpi_integer,
c & iset,1,mpi_integer,king,
#ifdef MPI
c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
c & " nfgtasks",nfgtasks
+ call flush(iout)
if (nfgtasks.gt.1) then
#ifdef MPI
time00=MPI_Wtime()
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
c call chainbuild_cart
endif
-c print *,'Processor',myrank,' calling etotal ipot=',ipot
+c write(iout,*) 'Processor',myrank,' calling etotal ipot=',ipot
c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
#else
c if (modecalc.eq.12.or.modecalc.eq.14) then
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+C BARTEK for dfa test!
+ if (wdfa_dist.gt.0) then
+ call edfad(edfadis)
+ else
+ edfadis=0
+ endif
+c print*, 'edfad is finished!', edfadis
+ if (wdfa_tor.gt.0) then
+ call edfat(edfator)
+ else
+ edfator=0
+ endif
+c print*, 'edfat is finished!', edfator
+ if (wdfa_nei.gt.0) then
+ call edfan(edfanei)
+ else
+ edfanei=0
+ endif
+c print*, 'edfan is finished!', edfanei
+ if (wdfa_beta.gt.0) then
+ call edfab(edfabet)
+ else
+ edfabet=0
+ endif
+c print*, 'edfab is finished!', edfabet
cmc
cmc Sep-06: egb takes care of dynamic ss bonds too
cmc
etors=0
edihcnstr=0
endif
+
+ if (constr_homology.ge.1) then
+ call e_modeller(ehomology_constr)
+ print *,'iset=',iset,'me=',me,ehomology_constr,
+ & 'Processor',fg_rank,' CG group',kolor,
+ & ' absolute rank',MyRank
+ else
+ ehomology_constr=0.0d0
+ endif
+
+
+c write(iout,*) ehomology_constr
c print *,"Processor",myrank," computed Utor"
C
C 6/23/01 Calculate double-torsional energy
C If performing constraint dynamics, call the constraint energy
C after the equilibration time
if(usampl.and.totT.gt.eq_time) then
+c write (iout,*) "CALL TO ECONSTR_BACK"
call EconstrQ
call Econstr_back
else
energia(21)=esccor
energia(22)=evdw_p
energia(23)=evdw_m
+ energia(24)=ehomology_constr
+ energia(25)=edfadis
+ energia(26)=edfator
+ energia(27)=edfanei
+ energia(28)=edfabet
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
if (dyn_ss) call dyn_set_nss
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ ehomology_constr=energia(24)
+ edfadis=energia(25)
+ edfator=energia(26)
+ edfanei=energia(27)
+ edfabet=energia(28)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+ & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+ & +wdfa_beta*edfabet
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+ & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+ & +wdfa_beta*edfabet
#endif
energia(0)=etot
c detecting NaNQ
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
enddo
enddo
#else
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
enddo
enddo
#endif
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
enddo
enddo
#endif
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ ehomology_constr=energia(24)
+C Bartek
+ edfadis = energia(25)
+ edfator = energia(26)
+ edfanei = energia(27)
+ edfabet = energia(28)
+
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
- & edihcnstr,ebr*nss,
- & Uconst,etot
+ & edihcnstr,ehomology_constr, ebr*nss,
+ & Uconst,edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei,
+ & edfabet,wdfa_beta,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
+ & 'EDFAD= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA distance energy)'/
+ & 'EDFAT= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA torsion energy)'/
+ & 'EDFAN= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA NCa energy)'/
+ & 'EDFAB= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA Beta energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ebr*nss,Uconst,etot
+ & ehomology_constr,ebr*nss,Uconst,edfadis,wdfa_dist,edfator,
+ & wdfa_tor,edfanei,wdfa_nei,edfabet,wdfa_beta,
+ & etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
+ & 'EDFAD= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA distance energy)'/
+ & 'EDFAT= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA torsion energy)'/
+ & 'EDFAN= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA NCa energy)'/
+ & 'EDFAB= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA Beta energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
& 'ebend',i,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
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)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
enddo
C Ufff.... We've done all this!!!
return
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
+ if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+ &(itype(i).eq.ntyp1)) cycle
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3) then
+ C if (i.gt.3) then
+ 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
+ ityp1=ithetyp(itype(i-2))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres) 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
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
common /sccalc/ time11,time12,time112,theti,it,nlobit
delta=0.02d0*pi
escloc=0.0D0
+c write(iout,*) "ESC: loc_start",loc_start," loc_end",loc_end
do i=loc_start,loc_end
costtab(i+1) =dcos(theta(i+1))
sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
return
end
c------------------------------------------------------------------------------
+c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA
+ subroutine e_modeller(ehomology_constr)
+ ehomology_constr=0.0
+ write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!"
+ return
+ end
+C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
+
+c------------------------------------------------------------------------------
subroutine etor_d(etors_d)
etors_d=0.0d0
return
return
end
c----------------------------------------------------------------------------
+c MODELLER restraint function
+ subroutine e_modeller(ehomology_constr)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+
+ integer nnn, i, j, k, ki, irec, l
+ integer katy, odleglosci, test7
+ real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
+ real*8 Eval,Erot
+ real*8 distance(max_template),distancek(max_template),
+ & min_odl,godl(max_template),dih_diff(max_template)
+
+c
+c FP - 30/10/2014 Temporary specifications for homology restraints
+c
+ double precision utheta_i,gutheta_i,sum_gtheta,sum_sgtheta,
+ & sgtheta
+ double precision, dimension (maxres) :: guscdiff,usc_diff
+ double precision, dimension (max_template) ::
+ & gtheta,dscdiff,uscdiffk,guscdiff2,guscdiff3,
+ & theta_diff
+c
+
+ 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'
+c
+c From subroutine Econstr_back
+c
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+c
+
+
+ do i=1,19
+ distancek(i)=9999999.9
+ enddo
+
+
+ odleg=0.0d0
+
+c Pseudo-energy and gradient from homology restraints (MODELLER-like
+c function)
+C AL 5/2/14 - Introduce list of restraints
+c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs start -------"
+#endif
+ do ii = link_start_homo,link_end_homo
+ i = ires_homo(ii)
+ j = jres_homo(ii)
+ dij=dist(i,j)
+c write (iout,*) "dij(",i,j,") =",dij
+ do k=1,constr_homology
+ distance(k)=odl(k,ii)-dij
+c write (iout,*) "distance(",k,") =",distance(k)
+ distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument
+c write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii)
+c write (iout,*) "distancek(",k,") =",distancek(k)
+c distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii)
+ enddo
+
+ min_odl=minval(distancek)
+c write (iout,* )"min_odl",min_odl
+#ifdef DEBUG
+ write (iout,*) "ij dij",i,j,dij
+ write (iout,*) "distance",(distance(k),k=1,constr_homology)
+ write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
+ write (iout,* )"min_odl",min_odl
+#endif
+ odleg2=0.0d0
+ do k=1,constr_homology
+c Nie wiem po co to liczycie jeszcze raz!
+c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
+c & (2*(sigma_odl(i,j,k))**2))
+ godl(k)=dexp(-distancek(k)+min_odl)
+ odleg2=odleg2+godl(k)
+
+ccc write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3,
+ccc & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=",
+ccc & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1),
+ccc & "sigma_odl(i,j,k)=", sigma_odl(i,j,k)
+
+ enddo
+c write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+c write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+ write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#endif
+ odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+c write (iout,*) "odleg",odleg ! sum of -ln-s
+c Gradient
+ sum_godl=odleg2
+ sum_sgodl=0.0d0
+ do k=1,constr_homology
+c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+c & *waga_dist)+min_odl
+c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
+ sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd
+ sum_sgodl=sum_sgodl+sgodl
+
+c sgodl2=sgodl2+sgodl
+c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1"
+c write(iout,*) "constr_homology=",constr_homology
+c write(iout,*) i, j, k, "TEST K"
+ enddo
+
+ if (homol_nset.gt.1)then
+ grad_odl3=waga_dist1(iset)*sum_sgodl/(sum_godl*dij)
+ else
+ grad_odl3=waga_dist*sum_sgodl/(sum_godl*dij)
+ endif
+c grad_odl3=sum_sgodl/(sum_godl*dij)
+
+
+c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2"
+c write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2),
+c & (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+
+ccc write(iout,*) godl, sgodl, grad_odl3
+
+c grad_odl=grad_odl+grad_odl3
+
+ do jik=1,3
+ ggodl=grad_odl3*(c(jik,i)-c(jik,j))
+ccc write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1))
+ccc write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl,
+ccc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+ ghpbc(jik,i)=ghpbc(jik,i)+ggodl
+ ghpbc(jik,j)=ghpbc(jik,j)-ggodl
+ccc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl,
+ccc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+c if (i.eq.25.and.j.eq.27) then
+c write(iout,*) "jik",jik,"i",i,"j",j
+c write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl
+c write(iout,*) "grad_odl3",grad_odl3
+c write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j)
+c write(iout,*) "ggodl",ggodl
+c write(iout,*) "ghpbc(",jik,i,")",
+c & ghpbc(jik,i),"ghpbc(",jik,j,")",
+c & ghpbc(jik,j)
+c endif
+ enddo
+ccc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=",
+ccc & dLOG(odleg2),"-odleg=", -odleg
+
+ enddo ! ii-loop for dist
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs end -------"
+c if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or.
+c & waga_d.eq.1.0d0) call sum_gradient
+#endif
+c Pseudo-energy and gradient from dihedral-angle restraints from
+c homology templates
+c write (iout,*) "End of distance loop"
+c call flush(iout)
+ kat=0.0d0
+c write (iout,*) idihconstr_start_homo,idihconstr_end_homo
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs start -------"
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ write (iout,*) "gloc_init(",i,icg,")",gloc(i,icg)
+ enddo
+#endif
+ 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)
+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)
+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)
+
+ kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
+ gdih(k)=dexp(kat3)
+ kat2=kat2+gdih(k)
+c write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
+c write(*,*)""
+ enddo
+c write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps
+c write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "i",i," betai",betai," kat2",kat2
+ write (iout,*) "gdih",(gdih(k),k=1,constr_homology)
+#endif
+ if (kat2.le.1.0d-14) cycle
+ kat=kat-dLOG(kat2/constr_homology)
+c write (iout,*) "kat",kat ! sum of -ln-s
+
+ccc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=",
+ccc & dLOG(kat2), "-kat=", -kat
+
+c ----------------------------------------------------------------------
+c Gradient
+c ----------------------------------------------------------------------
+
+ sum_gdih=kat2
+ sum_sgdih=0.0d0
+ do k=1,constr_homology
+ sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
+ sum_sgdih=sum_sgdih+sgdih
+ enddo
+c grad_dih3=sum_sgdih/sum_gdih
+ if (homol_nset.gt.1)then
+ grad_dih3=waga_angle1(iset)*sum_sgdih/sum_gdih
+ else
+ grad_dih3=waga_angle*sum_sgdih/sum_gdih
+ endif
+
+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
+c if (i.eq.25) then
+c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
+c endif
+ccc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3,
+ccc & gloc(nphi+i-3,icg)
+
+ enddo ! i-loop for dih
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs end -------"
+#endif
+
+c Pseudo-energy and gradient for theta angle restraints from
+c homology templates
+c FP 01/15 - inserted from econstr_local_test.F, loop structure
+c adapted
+
+c
+c For constr_homology reference structures (FP)
+c
+c Uconst_back_tot=0.0d0
+ Eval=0.0d0
+ Erot=0.0d0
+c Econstr_back legacy
+ do i=1,nres
+c do i=ithet_start,ithet_end
+ dutheta(i)=0.0d0
+c enddo
+c do i=loc_start,loc_end
+ do j=1,3
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
+ enddo
+c
+c do iref=1,nref
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "waga_theta",waga_theta
+ if (waga_theta.gt.0.0d0) then
+#ifdef DEBUG
+ write (iout,*) "usampl",usampl
+ write(iout,*) "------- theta restrs start -------"
+c do i=ithet_start,ithet_end
+c write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg)
+c enddo
+#endif
+c write (iout,*) "maxres",maxres,"nres",nres
+
+ do i=ithet_start,ithet_end
+c
+c do i=1,nfrag_back
+c ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset)
+c
+c Deviation of theta angles wrt constr_homology ref structures
+c
+ utheta_i=0.0d0 ! argument of Gaussian for single k
+ gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+c do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop
+c over residues in a fragment
+c write (iout,*) "theta(",i,")=",theta(i)
+ do k=1,constr_homology
+c
+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)
+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)
+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
+ enddo
+c write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps
+c write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps
+
+c
+c Gradient for multiple Gaussian restraint
+ sum_gtheta=gutheta_i
+ sum_sgtheta=0.0d0
+ do k=1,constr_homology
+c New generalized expr for multiple Gaussian from Econstr_back
+ sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd
+c
+c sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form?
+ sum_sgtheta=sum_sgtheta+sgtheta ! cum variable
+ enddo
+c grad_theta3=sum_sgtheta/sum_gtheta 1/*theta(i)? s. line below
+c grad_theta3=sum_sgtheta/sum_gtheta
+c
+c Final value of gradient using same var as in Econstr_back
+ dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
+c dutheta(i)=sum_sgtheta/sum_gtheta
+c
+c Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight
+ Eval=Eval-dLOG(gutheta_i/constr_homology)
+c write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps
+c write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s
+c Uconst_back=Uconst_back+utheta(i)
+ enddo ! (i-loop for theta)
+#ifdef DEBUG
+ write(iout,*) "------- theta restrs end -------"
+#endif
+ endif
+c
+c Deviation of local SC geometry
+c
+c Separation of two i-loops (instructed by AL - 11/3/2014)
+c
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c write (iout,*) "waga_d",waga_d
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs start -------"
+ write (iout,*) "Initial duscdiff,duscdiffx"
+ do i=loc_start,loc_end
+ write (iout,*) i,(duscdiff(jik,i),jik=1,3),
+ & (duscdiffx(jik,i),jik=1,3)
+ enddo
+#endif
+ do i=loc_start,loc_end
+ usc_diff_i=0.0d0 ! argument of Gaussian for single k
+ guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+c do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
+c write(iout,*) "xxtab, yytab, zztab"
+c write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)
+ do k=1,constr_homology
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ 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
+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,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
+c & xxref(j),yyref(j),zzref(j)
+ enddo
+c
+c Gradient
+c
+c Generalized expression for multiple Gaussian acc to that for a single
+c Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014)
+c
+c Original implementation
+c sum_guscdiff=guscdiff(i)
+c
+c sum_sguscdiff=0.0d0
+c do k=1,constr_homology
+c sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d?
+c sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff
+c sum_sguscdiff=sum_sguscdiff+sguscdiff
+c enddo
+c
+c Implementation of new expressions for gradient (Jan. 2015)
+c
+c grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !?
+ do k=1,constr_homology
+c
+c New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong
+c before. Now the drivatives should be correct
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+c
+c New implementation
+c
+ sum_guscdiff=guscdiff2(k)*!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong!
+ & sigma_d(k,i) ! for the grad wrt r'
+c sum_sguscdiff=sum_sguscdiff+sum_guscdiff
+c
+c
+c New implementation
+ sum_guscdiff = waga_d*sum_guscdiff
+ do jik=1,3
+ duscdiff(jik,i-1)=duscdiff(jik,i-1)+
+ & sum_guscdiff*(dXX_C1tab(jik,i)*dxx+
+ & dYY_C1tab(jik,i)*dyy+dZZ_C1tab(jik,i)*dzz)/guscdiff(i)
+ duscdiff(jik,i)=duscdiff(jik,i)+
+ & sum_guscdiff*(dXX_Ctab(jik,i)*dxx+
+ & dYY_Ctab(jik,i)*dyy+dZZ_Ctab(jik,i)*dzz)/guscdiff(i)
+ duscdiffx(jik,i)=duscdiffx(jik,i)+
+ & sum_guscdiff*(dXX_XYZtab(jik,i)*dxx+
+ & dYY_XYZtab(jik,i)*dyy+dZZ_XYZtab(jik,i)*dzz)/guscdiff(i)
+c
+#ifdef DEBUG
+ write(iout,*) "jik",jik,"i",i
+ write(iout,*) "dxx, dyy, dzz"
+ write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+ write(iout,*) "guscdiff2(",k,")",guscdiff2(k)
+c write(iout,*) "sum_sguscdiff",sum_sguscdiff
+cc write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i)
+c write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i)
+c write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i)
+c write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i)
+c write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i)
+c write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i)
+c write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i)
+c write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i)
+c write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i)
+c write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1)
+c write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i)
+c write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i)
+c endif
+#endif
+ enddo
+ enddo
+c
+c uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required?
+c usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ?
+c
+c write (iout,*) i," uscdiff",uscdiff(i)
+c
+c Put together deviations from local geometry
+
+c Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+
+c & wfrag_back(3,i,iset)*uscdiff(i)
+ Erot=Erot-dLOG(guscdiff(i)/constr_homology)
+c write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps
+c write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s
+c Uconst_back=Uconst_back+usc_diff(i)
+c
+c Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?)
+c
+c New implment: multiplied by sum_sguscdiff
+c
+
+ enddo ! (i-loop for dscdiff)
+
+c endif
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs end -------"
+ write (iout,*) "------ After SC loop in e_modeller ------"
+ do i=loc_start,loc_end
+ write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
+ write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
+ enddo
+ if (waga_theta.eq.1.0d0) then
+ write (iout,*) "in e_modeller after SC restr end: dutheta"
+ do i=ithet_start,ithet_end
+ write (iout,*) i,dutheta(i)
+ enddo
+ endif
+ if (waga_d.eq.1.0d0) then
+ write (iout,*) "e_modeller after SC loop: duscdiff/x"
+ do i=1,nres
+ write (iout,*) i,(duscdiff(j,i),j=1,3)
+ write (iout,*) i,(duscdiffx(j,i),j=1,3)
+ enddo
+ endif
+#endif
+
+c Total energy from homology restraints
+#ifdef DEBUG
+ write (iout,*) "odleg",odleg," kat",kat
+#endif
+c
+c Addition of energy of theta angle and SC local geom over constr_homologs ref strs
+c
+c ehomology_constr=odleg+kat
+ if (homol_nset.gt.1)then
+ ehomology_constr=waga_dist1(iset)*odleg+waga_angle1(iset)*kat+waga_theta*Eval
+ & +waga_d*Erot
+ else
+ ehomology_constr=waga_dist*odleg+waga_angle*kat+waga_theta*Eval
+ & +waga_d*Erot
+ endif
+c write (iout,*) "odleg",odleg," kat",kat," Uconst_back",Uconst_back
+c write (iout,*) "ehomology_constr",ehomology_constr
+c ehomology_constr=odleg+kat+Uconst_back
+ return
+c
+c FP 01/15 end
+c
+ 748 format(a8,f12.3,a6,f12.3,a7,f12.3)
+ 747 format(a12,i4,i4,i4,f8.3,f8.3)
+ 746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3)
+ 778 format(a7,1X,f10.3,1X,a4,1X,f10.3,1X,a5,1X,f10.3)
+ 779 format(i3,1X,i3,1X,i2,1X,a7,1X,f7.3,1X,a7,1X,f7.3,1X,a13,1X,
+ & f7.3,1X,a17,1X,f9.3,1X,a10,1X,f8.3,1X,a10,1X,f8.3)
+ end
+
+c------------------------------------------------------------------------------
subroutine etor_d(etors_d)
C 6/23/01 Compute double torsional energy
implicit real*8 (a-h,o-z)
esccor=0.0D0
do i=itau_start,itau_end
esccor_ii=0.0D0
+ if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
phii=phi(i)
projects / unres.git / commitdiff