From: Cezary Czaplewski Date: Sun, 22 Nov 2015 10:43:36 +0000 (+0100) Subject: CONSTR_HOMOL in mutlichain, debug and output cleaning X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=d0d389b657069f2239991163a15a46c42a36e097;hp=-c;p=unres.git CONSTR_HOMOL in mutlichain, debug and output cleaning --- d0d389b657069f2239991163a15a46c42a36e097 diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index b8a775e..0a868ce 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -1,14 +1,20 @@ 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 diff --git a/source/unres/src_MD-M/COMMON.MD b/source/unres/src_MD-M/COMMON.MD index c89ff32..131ea7c 100644 --- a/source/unres/src_MD-M/COMMON.MD +++ b/source/unres/src_MD-M/COMMON.MD @@ -13,6 +13,22 @@ & 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, @@ -28,14 +44,26 @@ 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 diff --git a/source/unres/src_MD-M/COMMON.REMD b/source/unres/src_MD-M/COMMON.REMD index 7109548..120731c 100644 --- a/source/unres/src_MD-M/COMMON.REMD +++ b/source/unres/src_MD-M/COMMON.REMD @@ -8,9 +8,10 @@ & 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, diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index da975d3..9ac2264 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -139,3 +139,6 @@ C Maximum number of conformation stored in cache on each CPU before sending 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) diff --git a/source/unres/src_MD-M/MREMD.F b/source/unres/src_MD-M/MREMD.F index 05d54ab..b3db8ab 100644 --- a/source/unres/src_MD-M/MREMD.F +++ b/source/unres/src_MD-M/MREMD.F @@ -1,3 +1,4 @@ +#define DEBUG subroutine MREMD implicit real*8 (a-h,o-z) include 'DIMENSIONS' @@ -29,7 +30,7 @@ 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 @@ -42,6 +43,7 @@ cold integer nup(0:maxprocs),ndown(0:maxprocs) 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 @@ -58,9 +60,19 @@ cdeb imin_itime_old=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) @@ -81,8 +93,9 @@ cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath 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 @@ -177,7 +190,7 @@ cd write (*,*) me," After broadcast: file_exist",file_exist 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,*) @@ -268,8 +281,13 @@ cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep) 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 @@ -551,7 +569,7 @@ C call enerprint(remd_ene(0,i)) 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" @@ -578,7 +596,7 @@ C call enerprint(remd_ene(0,i)) 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) @@ -616,6 +634,7 @@ c REMD - exchange & .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) @@ -625,7 +644,7 @@ cd call flush(iout) 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, @@ -693,21 +712,68 @@ c call mpi_gather(t_bath,1,mpi_double_precision, 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, @@ -753,10 +819,18 @@ cd end 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 @@ -770,20 +844,27 @@ co write(iout,*) 'REMD exchange temp,ene,elow,ehigh' 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 @@ -826,14 +907,14 @@ c write (iout,*) "ene_i_i",remd_ene(0,i) 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) @@ -844,7 +925,7 @@ c write (iout,*) "0,iex",remd_t_bath(iex) 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 @@ -949,9 +1030,9 @@ c call flush(iout) 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)) @@ -959,10 +1040,14 @@ cd write(iout,*) "########",ii 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 @@ -979,10 +1064,11 @@ cd write(iout,*) "i_dir=",i_dir 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 @@ -991,17 +1077,17 @@ cd write(iout,*) "i_dir=",i_dir 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. @@ -1013,33 +1099,39 @@ co & remd_t_bath(i) 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)) @@ -1054,7 +1146,7 @@ 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) @@ -1090,8 +1182,8 @@ cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx 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 @@ -1119,7 +1211,7 @@ c------------------------------------- & ,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, @@ -1147,10 +1239,17 @@ cd call flush(iout) & 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 @@ -1177,8 +1276,8 @@ co & " rescaling weights with temperature",t_bath 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 @@ -1231,6 +1330,7 @@ c----------------------------------------------------------------------- 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) @@ -1309,7 +1409,7 @@ c----------------------------------------------------------------------- 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) @@ -1370,7 +1470,7 @@ c----------------------------------------------------------------------- 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) @@ -1503,8 +1603,8 @@ c end debugging 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) @@ -1631,6 +1731,7 @@ c end debugging 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 @@ -1718,6 +1819,14 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) #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, @@ -1752,7 +1861,7 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) enddo - if(usampl) then + if(usampl.or.homol_nset.gt.1) then #ifdef AIX if(me.eq.king)then call xdrfint_(ixdrf, nset, iret) @@ -1794,7 +1903,7 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) 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, @@ -1802,6 +1911,11 @@ c & CG_COMM,ierr) 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 diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 91353c1..79cd1b9 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -214,6 +214,18 @@ cd print *,'nterm=',nterm 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 @@ -321,6 +333,7 @@ C 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" @@ -414,20 +427,23 @@ cMS$ATTRIBUTES C :: proc_proc 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 @@ -480,6 +496,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.TIME1' include 'COMMON.MAXGRAD' include 'COMMON.SCCOR' + include 'COMMON.MD' #ifdef TIMING time01=MPI_Wtime() #endif @@ -720,6 +737,14 @@ c enddo & +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 @@ -1006,6 +1031,7 @@ C------------------------------------------------------------------------ estr=energia(17) Uconst=energia(20) esccor=energia(21) + ehomology_constr=energia(24) eliptran=energia(22) Eafmforce=energia(23) #ifdef SPLITELE @@ -1015,7 +1041,7 @@ C------------------------------------------------------------------------ & 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)'/ @@ -1027,7 +1053,7 @@ C------------------------------------------------------------------------ & '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)'/ @@ -1038,6 +1064,7 @@ C------------------------------------------------------------------------ & '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)'/ @@ -1051,7 +1078,9 @@ C------------------------------------------------------------------------ & 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)'/ @@ -1072,6 +1101,7 @@ C------------------------------------------------------------------------ & '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)'/ @@ -3649,7 +3679,7 @@ cd & xmedi,ymedi,zmedi,xj,yj,zj 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 @@ -4514,8 +4544,8 @@ c i+3 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 @@ -5041,8 +5071,8 @@ c if (sss.eq.0) print *,'czasem jest OK' 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 @@ -6661,6 +6691,15 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) 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 @@ -6775,6 +6814,596 @@ cd write (iout,*) 'edihcnstr',edihcnstr 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) diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F index 0b711bd..09cc511 100644 --- a/source/unres/src_MD-M/geomout.F +++ b/source/unres/src_MD-M/geomout.F @@ -215,7 +215,7 @@ c------------------------------------------------------------------------ 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), diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F index acd4472..7295d2e 100644 --- a/source/unres/src_MD-M/gradient_p.F +++ b/source/unres/src_MD-M/gradient_p.F @@ -389,6 +389,17 @@ C 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 diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index 7ee3e42..b71e9e6 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -1461,5 +1461,49 @@ 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 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 183a917..efc7d64 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -1005,8 +1005,6 @@ c Dtilde(2,2,i)=0.0d0 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 @@ -1122,7 +1120,7 @@ C do j=1,ntyp 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 diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F index 2026425..a504ff7 100644 --- a/source/unres/src_MD-M/readpdb.F +++ b/source/unres/src_MD-M/readpdb.F @@ -13,16 +13,26 @@ C geometry. 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 @@ -66,53 +76,91 @@ C End current chain 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 @@ -124,7 +172,7 @@ C Calculate dummy residue coordinates inside the "chain" of a multichain 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 @@ -132,14 +180,12 @@ 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 - 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 @@ -179,7 +225,7 @@ C Calculate the CM of the last side chain. do j=1,3 dc(j,ires)=sccor(j,iii) enddo - else + else call sccenter(ires,iii,sccor) endif nsup=nres @@ -205,7 +251,7 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue 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 @@ -238,8 +284,28 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue 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), @@ -247,7 +313,7 @@ C Calculate internal coordinates. 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) @@ -290,7 +356,7 @@ c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) 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 @@ -309,7 +375,7 @@ c enddo c enddiagnostic C makes copy of chains write (iout,*) "symetr", symetr - + if (symetr.gt.1) then call permut(symetr) nperm=1 @@ -355,7 +421,7 @@ c diag 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 @@ -377,7 +443,7 @@ c--------------------------------------------------------------------------- include 'DIMENSIONS' #ifdef MPI include "mpif.h" -#endif +#endif include 'COMMON.LOCAL' include 'COMMON.VAR' include 'COMMON.CHAIN' @@ -387,7 +453,8 @@ c--------------------------------------------------------------------------- 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 @@ -400,11 +467,11 @@ c--------------------------------------------------------------------------- & '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 @@ -444,6 +511,9 @@ c endif 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 @@ -477,7 +547,7 @@ c------------------------------------------------------------------------------- include 'DIMENSIONS' #ifdef MPI include "mpif.h" -#endif +#endif include 'COMMON.LOCAL' include 'COMMON.VAR' include 'COMMON.CHAIN' @@ -600,4 +670,365 @@ c print *,vbld(i+1),vbld(i+1+nres) 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 diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 782e453..7be445d 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -98,6 +98,7 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file 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 @@ -139,7 +140,7 @@ C Set up the time limit (caution! The time must be input in minutes!) 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) @@ -766,6 +767,17 @@ C 12/1/95 Added weight for the multi-body term WCORR 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, @@ -795,6 +807,8 @@ c print *,'Finished reading pdb data' & 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 @@ -968,7 +982,6 @@ c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup 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 @@ -985,6 +998,56 @@ c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup 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 @@ -1904,8 +1967,16 @@ C Get parameter filenames and open the parameter files. 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) @@ -1930,9 +2001,17 @@ c print *,"Processor",myrank," opened file IBOND" 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" @@ -1961,8 +2040,16 @@ C Get parameter filenames and open the parameter files. 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) @@ -1995,9 +2082,9 @@ C Get parameter filenames and open the parameter files. 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) @@ -2209,6 +2296,7 @@ c------------------------------------------------------------------------------- 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 @@ -2218,7 +2306,7 @@ c------------------------------------------------------------------------------- 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) @@ -2236,33 +2324,36 @@ c------------------------------------------------------------------------------- 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 @@ -2280,7 +2371,7 @@ C--------------------------------------------------------------------------- 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) @@ -2333,6 +2424,18 @@ c write (iout,*) "IPAIR" 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 @@ -2425,6 +2528,413 @@ c write (iout,*) "j",j," k",k 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 @@ -2475,6 +2985,9 @@ C Initialize random number generator 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 @@ -2510,9 +3023,11 @@ C 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) @@ -2529,7 +3044,7 @@ C 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 diff --git a/source/unres/src_MD-M/unres.F b/source/unres/src_MD-M/unres.F index 2a828f2..03576c4 100644 --- a/source/unres/src_MD-M/unres.F +++ b/source/unres/src_MD-M/unres.F @@ -52,7 +52,7 @@ c call memmon_print_usage() 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