integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
- & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide
+ & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
+ & constr_homology,homol_nset
+ real*8 waga_homology
+ real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut
logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
& sideadd,lsecondary,read_cart,unres_pdb,
& vdisulf,searchsc,lmuca,dccart,extconf,out1file,
- & gnorm_check,gradout,split_ene
+ & gnorm_check,gradout,split_ene,read2sigma
common /cntrl/ modecalc,iscode,indpdb,indback,indphi,iranconf,
& icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint,
& overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
& ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
& constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog,
- & selfguide
+ & selfguide,
+ & constr_homology,homol_nset,read2sigma
+ common /homol/ waga_homology(maxprocs/20),
+ & waga_dist, waga_angle, waga_theta, waga_d, dist_cut
C... minim = .true. means DO minimization.
C... energy_dec = .true. means print energy decomposition matrix
& d_a(3,0:MAXRES2),d_a_work(6*MAXRES),kinetic_force(MAXRES6),
& Gsqrp(MAXRES2,MAXRES2),Gsqrm(MAXRES2,MAXRES2),
& vtot(MAXRES2),Gvec(maxres2,maxres2),Geigen(maxres2)
+
+ real*8 odl(max_template,maxdim),sigma_odl(max_template,maxdim),
+ & dih(max_template,maxres),sigma_dih(max_template,maxres),
+ & sigma_odlir(max_template,maxdim)
+c
+c Specification of new variables used in subroutine e_modeller
+c modified by FP (Nov.,2014)
+ real*8 xxtpl(max_template,maxres),yytpl(max_template,maxres),
+ & zztpl(max_template,maxres),thetatpl(max_template,maxres),
+ & sigma_theta(max_template,maxres),
+ & sigma_d(max_template,maxres)
+c
+
+ integer ires_homo(maxdim),
+ & jres_homo(maxdim),idomain(max_template,maxres)
+
double precision v_ini,d_time,d_time0,t_bath,tau_bath,
& EK,potE,potEcomp(0:n_ene+4),totE,totT,amax,kinetic_T,dvmax,damax,
& edriftmax,
integer n_timestep,ntwx,ntwe,lang,count_reset_moment,
& count_reset_vel,reset_fricmat,nfrag,npair,nfrag_back,
& ifrag_back(3,maxfrag_back,maxprocs/20),ntime_split,ntime_split0,
- & maxtime_split
+ & maxtime_split,lim_odl,lim_dih,link_start_homo,link_end_homo,
+ & idihconstr_start_homo,idihconstr_end_homo
logical large,print_compon,tbf,rest,reset_moment,reset_vel,
- & surfarea,rattle,usampl,mdpdb,RESPA,preminim
+ & surfarea,rattle,usampl,mdpdb,RESPA,preminim,
+ & l_homo(max_template,maxdim)
integer igmult_start,igmult_end,my_ng_count,ng_start,ng_counts,
& nginv_start,nginv_counts,myginv_ng_count
common /back_constr/ uconst_back,utheta,ugamma,uscdiff,
& dutheta,dugamma,duscdiff,duscdiffx,
& wfrag_back,nfrag_back,ifrag_back
+ common /homrestr/ odl,dih,sigma_dih,sigma_odl,
+ & lim_odl,lim_dih,ires_homo,jres_homo,link_start_homo,
+ & link_end_homo,idihconstr_start_homo,idihconstr_end_homo,
+ & idomain,l_homo
+c
+c FP (30/10/2014,04/03/2015)
+c
+ common /homrestr_double/
+ & xxtpl,yytpl,zztpl,thetatpl,sigma_theta,sigma_d,sigma_odlir
+c
common /qmeas/ qfrag,qpair,qinfrag,qinpair,wfrag,wpair,eq_time,
& Ucdfrag,Ucdpair,dUdconst,dUdxconst,dqwol,dxqwol,Uconst,
& iset,mset,nset,usampl,ifrag,ipair,npair,nfrag
& ndowna(0:maxprocs/4,0:maxprocs)
real t_restart1(5,maxprocs)
integer iset_restart1(maxprocs)
+ logical t_exchange_only
common /remdcommon/ nrep,nstex,retmin,retmax,remd_t,remd_tlist,
& remd_mlist,remd_m,mremdsync,restart1file,
- & traj1file,i_sync_step
+ & traj1file,i_sync_step,t_exchange_only
common /remdrestart/ i2rep,i2set,ifirst,nupa,ndowna,t_restart1,
& iset_restart1
real totT_cache,EK_cache,potE_cache,t_bath_cache,Uconst_cache,
C to master; depends on nstex / ntwx ratio
integer max_cache_traj
parameter (max_cache_traj=10)
+C Maximum number of templates in homology-modeling restraints
+ integer max_template
+ parameter(max_template=19)
+#define DEBUG
subroutine MREMD
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
integer iremd_iset(maxprocs)
integer*2 i_index
& (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
- double precision remd_ene(0:n_ene+4,maxprocs)
+ double precision remd_ene(0:n_ene+4,maxprocs),t_bath_old,e_tmp
integer iremd_acc(maxprocs),iremd_tot(maxprocs)
integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
integer ilen,rstcount
integer icache_all(maxprocs)
integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
logical synflag,end_of_run,file_exist /.false./,ovrtim
+ real ene_tol /1.0e-5/
cdeb imin_itime_old=0
ntwx_cache=0
endif
mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
-cd print *,'MREMD',nodes
+cd print *,'MREMD',nodes,homol_nset
cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
+ 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)
enddo
if(me.eq.king.or..not.out1file) then
- write(iout,*) (i2rep(i),i=0,nodes-1)
- write(iout,*) (i2set(i),i=0,nodes-1)
+ write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1)
+ write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
+ write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)"
do il=1,nset
do il1=1,mset(il)
do i=1,nrep
read (irest2,*) ndowna(0,il),
& (ndowna(i,il),i=1,ndowna(0,il))
enddo
- if(usampl) 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) then
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
iset=i2set(me)
+c broadcast iset to slaves
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ endif
if(me.eq.king.or..not.out1file)
& write(iout,*) me,"iset=",iset,"t_bath=",t_bath
endif
write (irest1,*) ndowna(0,il),
& (ndowna(i,il),i=1,ndowna(0,il))
enddo
- if(usampl) 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) then
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
write (irest2,*) iset
endif
close(irest2)
& .or.end_of_run.and.me.eq.king )
& .and. .not. mremdsync ) then
synflag=.true.
+ time01_=MPI_WTIME()
do i=1,nodes-1
call mpi_isend(itime,1,MPI_INTEGER,i,101,
& CG_COMM, ireqi(i), ierr)
call mpi_waitall(nodes-1,ireqi,statusi,ierr)
call mpi_barrier(CG_COMM, ierr)
time01=MPI_WTIME()
- write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
+ write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
if (out1file.or.traj1file) then
cdeb call mpi_gather(itime,1,mpi_integer,
cdeb & itime_all,1,mpi_integer,king,
c & remd_t_bath,1,mpi_double_precision,king,
c & CG_COMM,ierr)
potEcomp(n_ene+1)=t_bath
- if (usampl) then
+ t_bath_old=t_bath
+ if (usampl.or.homol_nset.gt.1) then
potEcomp(n_ene+2)=iset
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
+c broadcast iset to slaves and reduce energy
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call e_modeller(e_tmp)
+c write(iout,*) "iset+1 before reduce",e_tmp
+ call MPI_Barrier(FG_COMM,IERR)
+ call MPI_Reduce(e_tmp,potEcomp(n_ene+3),1,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ else
+ call e_modeller(potEcomp(n_ene+3))
+ endif
+c write(iout,*) "iset+1",potEcomp(n_ene+3)
+ else
+ call EconstrQ
+ potEcomp(n_ene+3)=Uconst
+ endif
iset=i_set_temp
+c broadcast iset to slaves
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ endif
+ else
+ potEcomp(n_ene+3)=0.0
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
+c broadcast iset to slaves and reduce energy
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call e_modeller(e_tmp)
+c write(iout,*) "iset-1 before reduce",e_tmp
+ call MPI_Barrier(FG_COMM,IERR)
+ call MPI_Reduce(e_tmp,potEcomp(n_ene+4),1,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ else
+ call e_modeller(potEcomp(n_ene+4))
+ endif
+c write(iout,*) "iset-1",potEcomp(n_ene+4)
+ else
+ call EconstrQ
+ potEcomp(n_ene+4)=Uconst
+ endif
iset=i_set_temp
+c broadcast iset to slaves
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ endif
+ else
+ potEcomp(n_ene+4)=0.0
endif
endif
call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
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
co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
do i=1,nodes
write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
enddo
endif
+#endif
c-------------------------------------
- IF(.not.usampl) 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
- call flush(iout)
+ctime call flush(iout)
write (iout,*) "remd_m(1)",remd_m(1)
+#endif
do irr=1,remd_m(1)
i=ifirst(iran_num(1,remd_m(1)))
+#ifdef DEBUG
write (iout,*) "i",i
- call flush(iout)
+#endif
+ctime call flush(iout)
do ii=1,nodes-1
+#ifdef DEBUG
write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
+#endif
if(i.gt.0.and.nupa(0,i).gt.0) then
iex=i
c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
c call flush(iout)
c write (iout,*) "rescaling weights with temperature",
c & remd_t_bath(iex)
- if (real(ene_i_i).ne.real(remd_ene(0,i))) then
+ if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
write (iout,*) "ERROR: inconsistent energies:",i,
& ene_i_i,remd_ene(0,i)
endif
call rescale_weights(remd_t_bath(iex))
c write (iout,*) "0,i",remd_t_bath(iex)
- call enerprint(remd_ene(0,i))
+c call enerprint(remd_ene(0,i))
call sum_energy(remd_ene(0,i),.false.)
c write (iout,*) "ene_i_iex",remd_ene(0,i)
c call enerprint(remd_ene(0,iex))
call sum_energy(remd_ene(0,iex),.false.)
- if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
+ if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
write (iout,*) "ERROR: inconsistent energies:",iex,
& ene_iex_iex,remd_ene(0,iex)
endif
enddo
cd write (iout,*) "exchange completed"
cd call flush(iout)
- ELSE
+ 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)
-cd write(iout,*) "i_dir=",i_dir
+ endif
+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
- call flush(iout)
+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
c following temperature changes.
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)
- & call mpi_scatter(iremd_iset,1,mpi_integer,
+ if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
+ call mpi_scatter(iremd_iset,1,mpi_integer,
& iset,1,mpi_integer,king,
& CG_COMM,ierr)
+c 8/31/2015 Correction by AL: send new iset to slaves
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ endif
+
+ endif
time07=MPI_WTIME()
if (me.eq.king .or. .not. out1file) then
cde write(iout,*) 'REMD after',me,t_bath
time08=MPI_WTIME()
if (me.eq.king .or. .not. out1file) then
- write(iout,*) 'REMD exchange time=',time08-time00
- call flush(iout)
+ write(iout,*) 'REMD exchange time=',time08-time02
+ctime call flush(iout)
endif
endif
enddo
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
#endif
enddo
enddo
+#ifdef DEBUG
+ write (iout,*) "Conformation read",il
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)')
+ & i,(d_restart1(j,i+2*nres*il),j=1,3),
+ & (d_restart1(j,nres+i+2*nres*il),j=1,3)
+ enddo
+#endif
enddo
endif
call mpi_scatter(d_restart1,3*2*nres,mpi_real,
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
-Corrected AL 8/19/2014: each processor needs whole iset array not only its
+c Corrected AL 8/19/2014: each processor needs whole iset array not only its
c own element
c call mpi_scatter(i2set,1,mpi_integer,
c & iset,1,mpi_integer,king,
call mpi_bcast(i2set(0),nodes,mpi_integer,king,
& CG_COMM,ierr)
iset=i2set(me)
+c broadcast iset to slaves
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ endif
endif
etors=0
edihcnstr=0
endif
+
+ if (constr_homology.ge.1) then
+ call e_modeller(ehomology_constr)
+c print *,'iset=',iset,'me=',me,ehomology_constr,
+c & 'Processor',fg_rank,' CG group',kolor,
+c & ' 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
energia(21)=esccor
energia(22)=eliptran
energia(23)=Eafmforce
+ energia(24)=ehomology_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"
esccor=energia(21)
eliptran=energia(22)
Eafmforce=energia(23)
+ ehomology_constr=energia(24)
#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+wliptran*eliptran+Eafmforce
+ & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+ & +wliptran*eliptran+Eafmforce
#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+wliptran*eliptran
+ & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+ & +wliptran*eliptran
& +Eafmforce
#endif
energia(0)=etot
include 'COMMON.TIME1'
include 'COMMON.MAXGRAD'
include 'COMMON.SCCOR'
+ include 'COMMON.MD'
#ifdef TIMING
time01=MPI_Wtime()
#endif
& +wliptran*gliptranx(j,i)
enddo
enddo
+ if (constr_homology.gt.0) then
+ do i=1,nct
+ do j=1,3
+ gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i)
+ gradx(j,i,icg)=gradx(j,i,icg)+duscdiffx(j,i)
+ enddo
+ enddo
+ endif
#ifdef DEBUG
write (iout,*) "gloc before adding corr"
do i=1,4*nres
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ ehomology_constr=energia(24)
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,ebr*nss,
+ & edihcnstr,ehomology_constr, ebr*nss,
& Uconst,eliptran,wliptran,Eafmforce,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
& 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
& 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
- & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6,
+ & 'EHPB= ',1pE16.6,' WEIGHT=',1pD16.6,
& ' (SS bridges & dist. cnstr.)'/
& 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
& 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
& '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)'/
& '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,
- & ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot
+ & ehomology_constr,ebr*nss,Uconst,
+ & eliptran,wliptran,Eafmforc,
+ & 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)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
if (energy_dec) then
write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)')
&'evdw1',i,j,evdwij
- &,iteli,itelj,aaa,evdw1
+c &,iteli,itelj,aaa,evdw1
write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
endif
eello_turn4=eello_turn4-(s1+s2+s3)
c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
- & 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
+c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
+c & 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
cd & ' eello_turn4_num',8*eello_turn4_num
#ifdef NEWCORR
evdwij=e1+e2
evdw2=evdw2+evdwij*sss
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
- & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
- & bad(itypj,iteli)
+ & 'evdw2',i,j,evdwij
+c & ,iteli,itypj,fac,aad(itypj,iteli),bad(itypj,iteli)
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
C
return
end
c------------------------------------------------------------------------------
+c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA
+ subroutine e_modeller(ehomology_constr)
+ ehomology_constr=0.0d0
+ 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
+c write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii)
+ if(.not.l_homo(k,ii)) cycle
+ distance(k)=odl(k,ii)-dij
+c write (iout,*) "distance(",k,") =",distance(k)
+c
+c For Gaussian-type Urestr
+c
+ 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)
+c
+c For Lorentzian-type Urestr
+c
+ if (waga_dist.lt.0.0d0) then
+ sigma_odlir(k,ii)=dsqrt(1/sigma_odl(k,ii))
+ distancek(k)=distance(k)**2/(sigma_odlir(k,ii)*
+ & (distance(k)**2+sigma_odlir(k,ii)**2))
+ endif
+ 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))
+ if(.not.l_homo(k,ii)) cycle
+ if (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ godl(k)=dexp(-distancek(k)+min_odl)
+ odleg2=odleg2+godl(k)
+c
+c For Lorentzian-type Urestr
+c
+ else
+ odleg2=odleg2+distancek(k)
+ endif
+
+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
+ if (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+c
+c For Lorentzian-type Urestr
+c
+ else
+ odleg=odleg+odleg2/constr_homology
+ endif
+c
+c write (iout,*) "odleg",odleg ! sum of -ln-s
+c Gradient
+c
+c For Gaussian-type Urestr
+c
+ if (waga_dist.ge.0.0d0) 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
+c
+ if(.not.l_homo(k,ii)) cycle
+ if (waga_dist.ge.0.0d0) then
+c For Gaussian-type Urestr
+c
+ sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd
+c
+c For Lorentzian-type Urestr
+c
+ else
+ sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+
+ & sigma_odlir(k,ii)**2)**2)
+ endif
+ 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 (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ grad_odl3=waga_homology(iset)*waga_dist
+ & *sum_sgodl/(sum_godl*dij)
+c
+c For Lorentzian-type Urestr
+c
+ else
+c Original grad expr modified by analogy w Gaussian-type Urestr grad
+c grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl
+ grad_odl3=-waga_homology(iset)*waga_dist*
+ & sum_sgodl/(constr_homology*dij)
+ endif
+c
+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
+ grad_dih3=waga_homology(iset)*waga_angle*sum_sgdih/sum_gdih
+
+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 Final value of gradient using same var as in Econstr_back
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)
+ & +sum_sgtheta/sum_gtheta*waga_theta
+ & *waga_homology(iset)
+c dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
+c & *waga_homology(iset)
+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_homology(iset)*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
+c
+c For Lorentzian-type Urestr
+c
+
+ if (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ ehomology_constr=(waga_dist*odleg+waga_angle*kat+
+ & waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c write (iout,*) "ehomology_constr=",ehomology_constr
+ else
+c
+c For Lorentzian-type Urestr
+c
+ ehomology_constr=(-waga_dist*odleg+waga_angle*kat+
+ & waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c write (iout,*) "ehomology_constr=",ehomology_constr
+ endif
+#ifdef DEBUG
+ write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat,
+ & "Eval",waga_theta,eval,
+ & "Erot",waga_d,Erot
+ write (iout,*) "ehomology_constr",ehomology_constr
+#endif
+ 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)
include 'COMMON.TORSION'
write (iout,'(/a)') 'Geometry of the virtual chain.'
write (iout,'(7a)') ' Res ',' d',' Theta',
- & ' Phi',' Dsc',' Alpha',' Omega'
+ & ' Gamma',' Dsc',' Alpha',' Beta'
do i=1,nres
iti=itype(i)
write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),
enddo
enddo
enddo
+c
+c Initialize the gradients of local restraints
+c
+ do i=1,nres
+ dutheta(i)=0.0d0
+ dugamma(i)=0.0d0
+ do j=1,3
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
+ enddo
C
C Initialize the gradient of local energy terms.
C
link_start=1
link_end=nhpb
#endif
+c write(2,*)"hpb_partition: link_start=",nhpb," link_end=",link_end
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine homology_partition
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.MD'
+ include 'COMMON.INTERACT'
+cd write(iout,*)"homology_partition: lim_odl=",lim_odl,
+cd & " lim_dih",lim_dih
+#ifdef MPI
+ if (me.eq.king .or. .not. out1file) write (iout,*) "MPI"
+ call int_bounds(lim_odl,link_start_homo,link_end_homo)
+ call int_bounds(lim_dih-nnt+1,idihconstr_start_homo,
+ & idihconstr_end_homo)
+ idihconstr_start_homo=idihconstr_start_homo+nnt-1
+ idihconstr_end_homo=idihconstr_end_homo+nnt-1
+ if (me.eq.king .or. .not. out1file)
+ & write (iout,*) 'Processor',fg_rank,' CG group',kolor,
+ & ' absolute rank',MyRank,
+ & ' lim_odl',lim_odl,' link_start=',link_start_homo,
+ & ' link_end',link_end_homo,' lim_dih',lim_dih,
+ & ' idihconstr_start_homo',idihconstr_start_homo,
+ & ' idihconstr_end_homo',idihconstr_end_homo
+#else
+ write (iout,*) "Not MPI"
+ link_start_homo=1
+ link_end_homo=lim_odl
+ idihconstr_start_homo=nnt
+ idihconstr_end_homo=lim_dih
+ write (iout,*)
+ & ' lim_odl',lim_odl,' link_start=',link_start_homo,
+ & ' link_end',link_end_homo,' lim_dih',lim_dih,
+ & ' idihconstr_start_homo',idihconstr_start_homo,
+ & ' idihconstr_end_homo',idihconstr_end_homo
+#endif
return
end
EEold(2,2,-i)=-b(10,i)+b(11,i)
EEold(2,1,-i)=-b(12,i)+b(13,i)
EEold(1,2,-i)=-b(12,i)-b(13,i)
- write(iout,*) "TU DOCHODZE"
- print *,"JESTEM"
c ee(1,1,i)=1.0d0
c ee(2,2,i)=1.0d0
c ee(2,1,i)=0.0d0
C epslip(i,j)=epslip(i,j)+0.05d0
C enddo
enddo
- write(iout,*) epslip(1,1),"OK?"
+c write(iout,*) epslip(1,1),"OK?"
C For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do i=1,ntyp
include 'COMMON.CONTROL'
include 'COMMON.DISTFIT'
include 'COMMON.SETUP'
- character*3 seq,atom,res
- character*80 card
- dimension sccor(3,20)
+ integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+ & ishift_pdb
+ logical lprn /.false./,fail
double precision e1(3),e2(3),e3(3)
+ double precision dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
+ character*80 card
+ double precision sccor(3,20)
integer rescode,iterter(maxres)
- logical fail
do i=1,maxres
iterter(i)=0
enddo
ibeg=1
+ ishift1=0
+ ishift=0
+c write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
lsecondary=.false.
nhfrag=0
nbfrag=0
endif
C Fish out the ATOM cards.
if (index(card(1:4),'ATOM').gt.0) then
- read (card(14:16),'(a3)') atom
- if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(12:16),*) atom
+c write (iout,*) "! ",atom," !",ires
+c if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+c write (iout,*) "ires",ires,ires-ishift+ishift1,
+c & " ires_old",ires_old
+c write (iout,*) "ishift",ishift," ishift1",ishift1
+c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
C Calculate the CM of the preceding residue.
if (ibeg.eq.0) then
if (unres_pdb) then
do j=1,3
- dc(j,ires+nres)=sccor(j,iii)
+ dc(j,ires)=sccor(j,iii)
enddo
else
- call sccenter(ires,iii,sccor)
+ call sccenter(ires_old,iii,sccor)
endif
+ iii=0
endif
C Start new residue.
-c write (iout,'(a80)') card
- read (card(23:26),*) ires
- read (card(18:20),'(a3)') res
- if (ibeg.eq.1) then
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+ write (iout,*) "BEG ires",ires
ishift=ires-1
if (res.ne.'GLY' .and. res.ne. 'ACE') then
ishift=ishift-1
itype(1)=ntyp1
endif
+ ires=ires-ishift+ishift1
+ ires_old=ires
+c write (iout,*) "ishift",ishift," ires",ires,
+c & " ires_old",ires_old
c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
ibeg=0
else if (ibeg.eq.2) then
c Start a new chain
ishift=-ires_old+ires-1
-c write (iout,*) "New chain started",ires,ishift
- ibeg=0
+ ires=ires_old+1
+ write (iout,*) "New chain started",ires,ishift
+ ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
endif
- ires=ires-ishift
-c write (2,*) "ires",ires," ishift",ishift
- if (res.eq.'ACE') then
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
itype(ires)=10
else
itype(ires)=rescode(ires,res,0)
endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+c write (iout,*) "ires_old",ires_old," ires",ires
+c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c ishift1=ishift1+1
+c endif
+c write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+ & res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
- if(me.eq.king.or..not.out1file)
- & write (iout,'(2i3,2x,a,3f8.3)')
- & ires,itype(ires),res,(c(j,ires),j=1,3)
- iii=1
+c write (iout,*) "backbone ",atom
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)')
+ & ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
do j=1,3
sccor(j,iii)=c(j,ires)
enddo
- else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
- & atom.ne.'N ' .and. atom.ne.'C ') then
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+c write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
+ & atom.ne.'N' .and. atom.ne.'C' .and.
+ & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
+ & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c write (iout,*) "sidechain ",atom
iii=iii+1
read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
endif
C system
nres=ires
do i=2,nres-1
- write (iout,*) i,itype(i),itype(i+1)
+c write (iout,*) i,itype(i),itype(i+1)
if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
if (unres_pdb) then
C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- print *,i,'tu dochodze'
call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
if (fail) then
e2(1)=0.0d0
e2(2)=1.0d0
e2(3)=0.0d0
endif !fail
- print *,i,'a tu?'
do j=1,3
c(j,i)=c(j,i-1)-1.9d0*e2(j)
enddo
do j=1,3
dc(j,ires)=sccor(j,iii)
enddo
- else
+ else
call sccenter(ires,iii,sccor)
endif
nsup=nres
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
- endif
+ endif
endif
do i=2,nres-1
do j=1,3
enddo
endif
endif
+C Copy the coordinates to reference coordinates
+c do i=1,2*nres
+c do j=1,3
+c cref(j,i)=c(j,i)
+c enddo
+c enddo
+C Calculate internal coordinates.
+ if (lprn) then
+ write (iout,'(/a)')
+ & "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+ & (c(j,ires+nres),j=1,3)
+ enddo
+ endif
C Calculate internal coordinates.
if(me.eq.king.or..not.out1file)then
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
do ires=1,nres
write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
& ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
enddo
endif
call int_from_cart(.true.,.false.)
- call sc_loc_geom(.true.)
+ call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
chain_rep(j,lll,kkk)=c(j,i)
chain_rep(j,lll+nres,kkk)=c(j,i+nres)
endif
- enddo
+ enddo
enddo
write (iout,*) chain_length
if (chain_length.eq.0) chain_length=nres
c enddiagnostic
C makes copy of chains
write (iout,*) "symetr", symetr
-
+
if (symetr.gt.1) then
call permut(symetr)
nperm=1
1 ' ', 6X,'X',11X,'Y',11X,'Z',
& 10X,'X',11X,'Y',11X,'Z')
110 format (a,'(',i3,')',6f12.5)
-
+
enddo
cc enddiag
do j=1,nbfrag
include 'DIMENSIONS'
#ifdef MPI
include "mpif.h"
-#endif
+#endif
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.NAMES'
include 'COMMON.CONTROL'
include 'COMMON.SETUP'
- character*3 seq,atom,res
+ character*3 seq,res
+c character*5 atom
character*80 card
dimension sccor(3,20)
integer rescode
& 'Internal coordinates calculated from crystal structure.'
if (lside) then
write (iout,'(8a)') ' Res ',' dvb',' Theta',
- & ' Phi',' Dsc_id',' Dsc',' Alpha',
- & ' Omega'
+ & ' Gamma',' Dsc_id',' Dsc',' Alpha',
+ & ' Beta '
else
write (iout,'(4a)') ' Res ',' dvb',' Theta',
- & ' Phi'
+ & ' Gamma'
endif
endif
#ifdef MPI
enddo
iti=itype(i)
di=dist(i,nres+i)
+C 10/03/12 Adam: Correction for zero SC-SC bond length
+ if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
+ & di=dsc(itype(i))
vbld(i+nres)=di
if (itype(i).ne.10) then
vbld_inv(i+nres)=1.0d0/di
include 'DIMENSIONS'
#ifdef MPI
include "mpif.h"
-#endif
+#endif
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
enddo
return
end
-
+
+ subroutine readpdb_template(k)
+C Read the PDB file for read_constr_homology with read2sigma
+C and convert the peptide geometry into virtual-chain geometry.
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.LOCAL'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.NAMES'
+ include 'COMMON.CONTROL'
+ include 'COMMON.DISTFIT'
+ include 'COMMON.SETUP'
+ integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+ & ishift_pdb
+ logical lprn /.false./,fail
+ double precision e1(3),e2(3),e3(3)
+ double precision dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
+ character*80 card
+ double precision sccor(3,20)
+ integer rescode,iterter(maxres)
+ do i=1,maxres
+ iterter(i)=0
+ enddo
+ ibeg=1
+ ishift1=0
+ ishift=0
+c write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
+ do
+ read (ipdbin,'(a80)',end=10) card
+ if (card(:3).eq.'END') then
+ goto 10
+ else if (card(:3).eq.'TER') then
+C End current chain
+ ires_old=ires+2
+ itype(ires_old-1)=ntyp1
+ iterter(ires_old-1)=1
+ itype(ires_old)=ntyp1
+ iterter(ires_old)=1
+ ibeg=2
+ write (iout,*) "Chain ended",ires,ishift,ires_old
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+C Fish out the ATOM cards.
+ if (index(card(1:4),'ATOM').gt.0) then
+ read (card(12:16),*) atom
+c write (iout,*) "! ",atom," !",ires
+c if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+c write (iout,*) "ires",ires,ires-ishift+ishift1,
+c & " ires_old",ires_old
+c write (iout,*) "ishift",ishift," ishift1",ishift1
+c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
+C Calculate the CM of the preceding residue.
+ if (ibeg.eq.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires_old,iii,sccor)
+ endif
+ iii=0
+ endif
+C Start new residue.
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+c write (iout,*) "BEG ires",ires
+ ishift=ires-1
+ if (res.ne.'GLY' .and. res.ne. 'ACE') then
+ ishift=ishift-1
+ itype(1)=ntyp1
+ endif
+ ires=ires-ishift+ishift1
+ ires_old=ires
+c write (iout,*) "ishift",ishift," ires",ires,
+c & " ires_old",ires_old
+c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ibeg=0
+ else if (ibeg.eq.2) then
+c Start a new chain
+ ishift=-ires_old+ires-1
+ ires=ires_old+1
+ write (iout,*) "New chain started",ires,ishift
+ ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
+ endif
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
+ itype(ires)=10
+ else
+ itype(ires)=rescode(ires,res,0)
+ endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+c write (iout,*) "ires_old",ires_old," ires",ires
+c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c ishift1=ishift1+1
+c endif
+c write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+ & res.eq.'NHE'.and.atom(:2).eq.'HN') then
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)')
+ & ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
+ do j=1,3
+ sccor(j,iii)=c(j,ires)
+ enddo
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+c write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
+ & atom.ne.'N' .and. atom.ne.'C' .and.
+ & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
+ & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c write (iout,*) "sidechain ",atom
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+ endif
+ enddo
+ 10 if(me.eq.king.or..not.out1file)
+ & write (iout,'(a,i5)') ' Nres: ',ires
+C Calculate dummy residue coordinates inside the "chain" of a multichain
+C system
+ nres=ires
+ do i=2,nres-1
+c write (iout,*) i,itype(i),itype(i+1)
+ if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
+ if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
+C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif !fail
+ do j=1,3
+ c(j,i)=c(j,i-1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i-2)-c(j,i-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i-1)+dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ else !itype(i+1).eq.ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i+3)-c(j,i+2))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i+1)-dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ endif !itype(i+1).eq.ntyp1
+ endif !itype.eq.ntyp1
+ enddo
+C Calculate the CM of the last side chain.
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ nsup=nres
+ nstart_sup=1
+ if (itype(nres).ne.10) then
+ nres=nres+1
+ itype(nres)=ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,nres)=c(j,nres-1)+dcj
+ c(j,2*nres)=c(j,nres)
+ enddo
+ endif
+ endif
+ do i=2,nres-1
+ do j=1,3
+ c(j,i+nres)=dc(j,i)
+ enddo
+ enddo
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ enddo
+ if (itype(1).eq.ntyp1) then
+ nsup=nsup-1
+ nstart_sup=2
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,4)-c(j,3))/2.0
+ c(j,1)=c(j,2)-dcj
+ c(j,nres+1)=c(j,1)
+ enddo
+ endif
+ endif
+C Copy the coordinates to reference coordinates
+c do i=1,2*nres
+c do j=1,3
+c cref(j,i)=c(j,i)
+c enddo
+c enddo
+C Calculate internal coordinates.
+ if (lprn) then
+ write (iout,'(/a)')
+ & "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+ & (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+C Calculate internal coordinates.
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
+ & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+ & (c(j,nres+ires),j=1,3)
+ enddo
+ endif
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+c & vbld_inv(i+nres)
+ enddo
+c call chainbuild
+C Copy the coordinates to reference coordinates
+C Splits to single chain if occurs
+ kkk=1
+ lll=0
+ cou=1
+ do i=1,nres
+ lll=lll+1
+cc write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
+ if (i.gt.1) then
+ if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
+ chain_length=lll-1
+ kkk=kkk+1
+c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
+ lll=1
+ endif
+ endif
+ do j=1,3
+ cref(j,i,cou)=c(j,i)
+ cref(j,i+nres,cou)=c(j,i+nres)
+ if (i.le.nres) then
+ chain_rep(j,lll,kkk)=c(j,i)
+ chain_rep(j,lll+nres,kkk)=c(j,i+nres)
+ endif
+ enddo
+ enddo
+ write (iout,*) chain_length
+ if (chain_length.eq.0) chain_length=nres
+ do j=1,3
+ chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
+ chain_rep(j,chain_length+nres,symetr)
+ &=chain_rep(j,chain_length+nres,1)
+ enddo
+c diagnostic
+c write (iout,*) "spraw lancuchy",chain_length,symetr
+c do i=1,4
+c do kkk=1,chain_length
+c write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
+c enddo
+c enddo
+c enddiagnostic
+
+ return
+ end
C Set up the time limit (caution! The time must be input in minutes!)
read_cart=index(controlcard,'READ_CART').gt.0
call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
call readi(controlcard,'SYM',symetr,1)
call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
extconf=(index(controlcard,'EXTCONF').gt.0)
AFMlog=(index(controlcard,'AFM'))
selfguide=(index(controlcard,'SELFGUIDE'))
- print *,'AFMlog',AFMlog,selfguide,"KUPA"
+c print *,'AFMlog',AFMlog,selfguide,"KUPA"
call readi(controlcard,'IPRINT',iprint,0)
call readi(controlcard,'MAXGEN',maxgen,10000)
call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
34 continue
c print *,'Begin reading pdb data'
call readpdb
+ do i=1,2*nres
+ do j=1,3
+ crefjlee(j,i)=c(j,i)
+ enddo
+ enddo
+#ifdef DEBUG
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
+ & (crefjlee(j,i+nres),j=1,3)
+ enddo
+#endif
c print *,'Finished reading pdb data'
if(me.eq.king.or..not.out1file)
& write (iout,'(a,i3,a,i3)')'nsup=',nsup,
& i,' after ',nsi,' trials'
endif
enddo
+C 10/03/12 Adam: Recalculate coordinates with new side chain positions
+ call chainbuild
endif
endif
if (indpdb.eq.0) then
if (constr_dist.gt.0) call read_dist_constr
write (iout,*) "After read_dist_constr nhpb",nhpb
if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
- call hpb_partition
if(me.eq.king.or..not.out1file)
& write (iout,*) 'Contact order:',co
if (pdbref) then
enddo
endif
endif
+
+
+ if (constr_homology.gt.0) then
+ call read_constr_homology
+ if (indpdb.gt.0 .or. pdbref) then
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=crefjlee(j,i)
+ cref(j,i,1)=crefjlee(j,i)
+ enddo
+ enddo
+ endif
+#ifdef DEBUG
+ write (iout,*) "Array C"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
+ & (c(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Array Cref"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
+ & (cref(j,i+nres,1),j=1,3)
+ enddo
+#endif
+ call int_from_cart1(.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ enddo
+ else
+ homol_nset=0
+ endif
+
+
+ if (nhpb.gt.0) call hpb_partition
+c write (iout,*) "After read_dist_constr nhpb",nhpb
+c call flush(iout)
if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
& .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
& modecalc.ne.10) then
open (ibond,file=bondname,status='old',readonly,shared)
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old',readonly,shared)
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
+#endif
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',readonly,shared)
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',readonly,shared)
call getenv_loc('TORDPAR',tordname)
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old',action='read')
c print *,"Processor",myrank," opened file ITHEP"
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old',action='read')
+#endif
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
c print *,"Processor",myrank," opened file IROTAM"
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',action='read')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
c print *,"Processor",myrank," opened file ITORP"
open (ibond,file=bondname,status='old')
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old')
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old')
+#endif
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old')
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old')
call getenv_loc('TORDPAR',tordname)
open (isccor,file=sccorname,status='old',readonly)
#ifndef CRYST_THETA
call getenv_loc('THETPARPDB',thetname_pdb)
- print *,"thetname_pdb ",thetname_pdb
+c print *,"thetname_pdb ",thetname_pdb
open (ithep_pdb,file=thetname_pdb,status='old',action='read')
- print *,ithep_pdb," opened"
+c print *,ithep_pdb," opened"
#endif
call getenv_loc('FOURIER',fouriername)
open (ifourier,file=fouriername,status='old',readonly)
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.MD'
+ include 'COMMON.CONTROL'
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
totTafm=totT
do i=1,2*nres
read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
read (irest2,*) iset
endif
close(irest2)
include 'COMMON.IOUNITS'
include 'COMMON.MD'
include 'COMMON.CONTROL'
+ integer iset1
read(inp,*) nset,nfrag,npair,nfrag_back
if(me.eq.king.or..not.out1file)
& write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
& " nfrag_back",nfrag_back
- do iset=1,nset
- read(inp,*) mset(iset)
+ do iset1=1,nset
+ read(inp,*) mset(iset1)
do i=1,nfrag
- read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
- & qinfrag(i,iset)
+ read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
+ & qinfrag(i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
- & ifrag(2,i,iset), qinfrag(i,iset)
+ & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
+ & ifrag(2,i,iset1), qinfrag(i,iset1)
enddo
do i=1,npair
- read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
- & qinpair(i,iset)
+ read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
+ & qinpair(i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
- & ipair(2,i,iset), qinpair(i,iset)
+ & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
+ & ipair(2,i,iset1), qinpair(i,iset1)
enddo
do i=1,nfrag_back
- read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
- & wfrag_back(3,i,iset),
- & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
+ & wfrag_back(3,i,iset1),
+ & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
- & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ & write(iout,*) "A",i,wfrag_back(1,i,iset1),
+ & wfrag_back(2,i,iset1),
+ & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
+ & ifrag_back(2,i,iset1)
enddo
enddo
return
include 'COMMON.IOUNITS'
include 'COMMON.SBRIDGE'
character*320 afmcard
- print *, "wchodze"
+c print *, "wchodze"
call card_concat(afmcard)
call readi(afmcard,"BEG",afmbeg,0)
call readi(afmcard,"END",afmend,0)
c do i=1,npair_
c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
c enddo
+ if (.not.refstr .and. nfrag.gt.0) then
+ write (iout,*)
+ & "ERROR: no reference structure to compute distance restraints"
+ write (iout,*)
+ & "Restraints must be specified explicitly (NDIST=number)"
+ stop
+ endif
+ if (nfrag.lt.2 .and. npair.gt.0) then
+ write (iout,*) "ERROR: Less than 2 fragments specified",
+ & " but distance restraints between pairs requested"
+ stop
+ endif
call flush(iout)
do i=1,nfrag_
if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
return
end
c-------------------------------------------------------------------------------
+
+ subroutine read_constr_homology
+
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.GEO'
+ include 'COMMON.INTERACT'
+c
+c For new homol impl
+c
+ include 'COMMON.VAR'
+c
+
+c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
+c & dist_cut
+c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+c & sigma_odl_temp(maxres,maxres,max_template)
+ character*2 kic2
+ character*24 model_ki_dist, model_ki_angle
+ character*500 controlcard
+ integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
+ logical lprn /.true./
+ integer ilen
+ external ilen
+c
+c FP - Nov. 2014 Temporary specifications for new vars
+c
+ double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+ double precision, dimension (max_template,maxres) :: rescore
+ double precision, dimension (max_template,maxres) :: rescore2
+ character*24 pdbfile,tpl_k_rescore
+c -----------------------------------------------------------------
+c Reading multiple PDB ref structures and calculation of retraints
+c not using pre-computed ones stored in files model_ki_{dist,angle}
+c FP (Nov., 2014)
+c -----------------------------------------------------------------
+c
+c
+c Alternative: reading from input
+ call card_concat(controlcard)
+ call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
+ call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
+ call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
+ call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
+
+ call readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
+ if (homol_nset.gt.1)then
+ call card_concat(controlcard)
+ read(controlcard,*) (waga_homology(i),i=1,homol_nset)
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "iset homology_weight "
+ do i=1,homol_nset
+ write(iout,*) i,waga_homology(i)
+ enddo
+ endif
+ iset=mod(kolor,homol_nset)+1
+ else
+ iset=1
+ waga_homology(1)=1.0
+ endif
+
+cd write (iout,*) "nnt",nnt," nct",nct
+cd call flush(iout)
+
+
+ lim_odl=0
+ lim_dih=0
+c
+c New
+c
+ lim_theta=0
+ lim_xx=0
+c
+ write(iout,*) 'nnt=',nnt,'nct=',nct
+c
+ do i = nnt,nct
+ do k=1,constr_homology
+ idomain(k,i)=0
+ enddo
+ enddo
+
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ ii_in_use(ii)=0
+ enddo
+ enddo
+
+ do k=1,constr_homology
+
+ read(inp,'(a)') pdbfile
+c Next stament causes error upon compilation (?)
+c if(me.eq.king.or. .not. out1file)
+c write (iout,'(2a)') 'PDB data will be read from file ',
+c & pdbfile(:ilen(pdbfile))
+ write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ stop
+ 34 continue
+c print *,'Begin reading pdb data'
+c
+c Files containing res sim or local scores (former containing sigmas)
+c
+
+ write(kic2,'(bz,i2.2)') k
+
+ tpl_k_rescore="template"//kic2//".sco"
+
+ unres_pdb=.false.
+ if (read2sigma) then
+ call readpdb_template(k)
+ else
+ call readpdb
+ endif
+c
+c Distance restraints
+c
+c ... --> odl(k,ii)
+C Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=cref(j,i,1)
+c write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+c
+c From read_dist_constr (commented out 25/11/2014 <-> res sim)
+c
+c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
+ open (ientin,file=tpl_k_rescore,status='old')
+ if (nnt.gt.1) rescore(k,1)=0.0d0
+ do irec=nnt,maxdim ! loop for reading res sim
+ if (read2sigma) then
+ read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
+ & idomain_tmp
+ i_tmp=i_tmp+nnt-1
+ idomain(k,i_tmp)=idomain_tmp
+ rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0)
+ rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0)
+ else
+ idomain(k,irec)=1
+ read (ientin,*,end=1401) rescore_tmp
+
+c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
+ rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
+c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
+ endif
+ enddo
+ 1401 continue
+ close (ientin)
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+
+ if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
+
+ ii=ii+1
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+
+c write (iout,*) "k",k
+c write (iout,*) "i",i," j",j," constr_homology",
+c & constr_homology
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+c
+c Attempt to replace dist(i,j) by its definition in ...
+c
+ x12=c(1,i)-c(1,j)
+ y12=c(2,i)-c(2,j)
+ z12=c(3,i)-c(3,j)
+ distal=dsqrt(x12*x12+y12*y12+z12*z12)
+ odl(k,ii)=distal
+ if (read2sigma) then
+ sigma_odl(k,ii)=0
+ do ik=i,j
+ sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+ enddo
+ sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+ if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
+ & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+ else
+ if (odl(k,ii).le.dist_cut) then
+ sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
+c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+ else
+#ifdef OLDSIGMA
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+#else
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+#endif
+ endif
+ endif
+c Following expr replaced by a positive exp argument
+c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+
+c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
+c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+c
+ sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
+c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
+c
+c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
+c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
+ else
+ ii=ii+1
+ l_homo(k,ii)=.false.
+ endif
+ enddo
+ enddo
+ lim_odl=ii
+ endif
+c
+c Theta, dihedral and SC retraints
+c
+ if (waga_angle.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_dih,status='old')
+c do irec=1,maxres-3 ! loop for reading sigma_dih
+c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
+c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
+c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_dih(k,i+nnt-1)
+c enddo
+c1402 continue
+c close (ientin)
+ do i = nnt+3,nct
+ if (idomain(k,i).eq.0) then
+ sigma_dih(k,i)=0.0
+ cycle
+ endif
+ dih(k,i)=phiref(i) ! right?
+c read (ientin,*) sigma_dih(k,i) ! original variant
+c write (iout,*) "dih(",k,i,") =",dih(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,i-1,") =",rescore(k,i-1),
+c & "rescore(",k,i-2,") =",rescore(k,i-2),
+c & "rescore(",k,i-3,") =",rescore(k,i-3)
+
+ sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)+rescore(k,i-3)
+ if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
+c write (iout,*) "Raw sigmas for dihedral angle restraints"
+c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
+c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
+c Instead of res sim other local measure of b/b str reliability possible
+ sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
+ enddo
+ lim_dih=nct-nnt-2
+ endif
+
+ if (waga_theta.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_theta,status='old')
+c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
+c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_theta(k,i+nnt-1)
+c enddo
+c1403 continue
+c close (ientin)
+
+ do i = nnt+2,nct ! right? without parallel.
+c do i = i=1,nres ! alternative for bounds acc to readpdb?
+c do i=ithet_start,ithet_end ! with FG parallel.
+ if (idomain(k,i).eq.0) then
+ sigma_theta(k,i)=0.0
+ cycle
+ endif
+ thetatpl(k,i)=thetaref(i)
+c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,i-1,") =",rescore(k,i-1),
+c & "rescore(",k,i-2,") =",rescore(k,i-2)
+c read (ientin,*) sigma_theta(k,i) ! 1st variant
+ sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)
+ if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
+ sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+
+c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+c rescore(k,i-2) ! right expression ?
+c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
+ enddo
+ endif
+ lim_theta=nct-nnt-1
+
+ if (waga_d.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_d,status='old')
+c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
+c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_d(k,i+nnt-1)
+c enddo
+c1404 continue
+
+ do i = nnt,nct ! right? without parallel.
+c do i=2,nres-1 ! alternative for bounds acc to readpdb?
+c do i=loc_start,loc_end ! with FG parallel.
+ if (itype(i).eq.10) cycle
+ if (idomain(k,i).eq.0 ) then
+ sigma_d(k,i)=0.0
+ cycle
+ endif
+ xxtpl(k,i)=xxref(i)
+ yytpl(k,i)=yyref(i)
+ zztpl(k,i)=zzref(i)
+c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
+c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
+c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
+ sigma_d(k,i)=rescore(k,i) ! right expression ?
+ sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+
+c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
+c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
+c read (ientin,*) sigma_d(k,i) ! 1st variant
+ if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
+ enddo
+ lim_xx=nct-nnt+1
+ endif
+ enddo
+c
+c remove distance restraints not used in any model from the list
+c shift data in all arrays
+c
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i=nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ if (ii_in_use(ii).eq.0) then
+ do ki=ii,lim_odl-1
+ ires_homo(ki)=ires_homo(ki+1)
+ jres_homo(ki)=jres_homo(ki+1)
+ ii_in_use(ki)=ii_in_use(ki+1)
+ do k=1,constr_homology
+ odl(k,ki)=odl(k,ki+1)
+ sigma_odl(k,ki)=sigma_odl(k,ki+1)
+ l_homo(k,ki)=l_homo(k,ki+1)
+ enddo
+ enddo
+ ii=ii-1
+ lim_odl=lim_odl-1
+ endif
+ enddo
+ enddo
+ endif
+ if (constr_homology.gt.0) call homology_partition
+ if (constr_homology.gt.0) call init_int_table
+cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
+cd & "lim_xx=",lim_xx
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c
+c Print restraints
+c
+ if (.not.lprn) return
+cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write (iout,*) "Distance restraints from templates"
+ do ii=1,lim_odl
+ write(iout,'(3i5,10(2f8.2,1x,l1,4x))')
+ & ii,ires_homo(ii),jres_homo(ii),
+ & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
+ & ki=1,constr_homology)
+ enddo
+ write (iout,*) "Dihedral angle restraints from templates"
+ do i=nnt+3,lim_dih
+ write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+ & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "Virtual-bond angle restraints from templates"
+ do i=nnt+2,lim_theta
+ write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+ & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "SC restraints from templates"
+ do i=nnt,lim_xx
+ write(iout,'(i5,10(4f8.2,4x))') i,
+ & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
+ & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
+ enddo
+ endif
+c -----------------------------------------------------------------
+ return
+ end
+c----------------------------------------------------------------------
+
#ifdef WINIFL
subroutine flush(iu)
return
C
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
+#ifdef AMD64
+ integer*8 iseedi8
+#endif
#ifdef MPI
include 'mpif.h'
logical OKRandom, prng_restart
if (fg_rank.eq.0) then
seed=seed*(me+1)+1
#ifdef AMD64
- if(me.eq.king)
- & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
- OKRandom = prng_restart(me,iseed)
+ iseedi8=dint(seed)
+ if(me.eq.king .or. .not. out1file)
+ & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
+ write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
+ OKRandom = prng_restart(me,iseedi8)
#else
do i=1,4
tmp=65536.0d0**(4-i)
if (OKRandom) then
c r1 = prng_next(me)
r1=ran_number(0.0D0,1.0D0)
- if(me.eq.king)
+ if(me.eq.king .or. .not. out1file)
& write (iout,*) 'ran_num',r1
if (r1.lt.0.0d0) OKRandom=.false.
endif
call init_task
if (me.eq.king)
- & write(iout,*)'### LAST MODIFIED 4/25/08 7:29PM by adam'
+ & write(iout,*)'### LAST MODIFIED 11/22/15 10:36PM by czarek'
if (me.eq.king) call cinfo
C Read force field parameters and job setup data
call readrtns
projects / unres.git / commitdiff