& .and. mod(itime,count_reset_vel).eq.0) then
call random_vel
write(iout,'(a,f20.2)')
- & "Velocities reset to random values, time",totT
+ & "Velocities reset to random values, time",totT
do i=0,2*nres
do j=1,3
d_t_old(j,i)=d_t(j,i)
enddo
enddo
endif
- if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
+ if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
call inertia_tensor
call vcm_vel(vcm)
do j=1,3
call kinetic(EK)
kinetic_T=2.0d0/(dimen3*Rb)*EK
scalfac=dsqrt(T_bath/kinetic_T)
- write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
+ write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
do i=0,2*nres
do j=1,3
d_t_old(j,i)=scalfac*d_t(j,i)
stop
#endif
endif
+ itime_mat=itime
if (ntwe.ne.0) then
if (mod(itime,ntwe).eq.0) then
call statout(itime)
if (rstcount.eq.1000.or.itime.eq.n_timestep) then
open(irest2,file=rest2name,status='unknown')
write(irest2,*) totT,EK,potE,totE,t_bath
- do i=1,2*nres
+ do i=0,2*nres-1
write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
enddo
- do i=1,2*nres
+ do i=0,2*nres-1
write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
close(irest2)
d_t_new(j,0)=d_t_old(j,0)+adt2
d_t(j,0)=d_t_old(j,0)+adt
enddo
- do i=nnt,nct-1
+ do i=nnt,nct-1
C SPYTAC ADAMA
C do i=0,nres
+#ifdef DEBUG
+ write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3)
+#endif
do j=1,3
adt=d_a_old(j,i)*d_time
adt2=0.5d0*adt
call fivediaginv_mult(dimen,fric_work, d_af_work)
c write (iout,*) "stochastic acceleratios"
call fivediaginv_mult(dimen,stochforcvec, d_as_work)
-c write (iout,*) "Leaving sddir_precalc"
#else
call ginv_mult(fric_work, d_af_work)
call ginv_mult(stochforcvec, d_as_work)
write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3)
write (iout,*) "d_as_work"
write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3)
+ write (iout,*) "Leaving sddir_precalc"
#endif
return
end
include 'COMMON.NAMES'
include 'COMMON.REMD'
include 'COMMON.TIME1'
+#ifdef LBFGS
+ character*9 status
+ integer niter
+ common /lbfgstat/ status,niter,nfun
+#endif
+ integer n_model_try,list_model_try(max_template),k
double precision tt0
real*8 energia_long(0:n_ene),
& energia_short(0:n_ene),vcm(3),incr(3)
integer iran_num
double precision etot
logical fail
+ integer i_start_models(0:nodes-1)
write (iout,*) "init_MD INDPDB",indpdb
d_time0=d_time
c write(iout,*) "d_time", d_time
if (me.eq.king)
& inquire(file=mremd_rst_name,exist=file_exist)
#ifdef MPI
- write (*,*) me," Before broadcast: file_exist",file_exist
+c write (*,*) me," Before broadcast: file_exist",file_exist
call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
& IERR)
- write (*,*) me," After broadcast: file_exist",file_exist
+c write (*,*) me," After broadcast: file_exist",file_exist
c inquire(file=mremd_rst_name,exist=file_exist)
#endif
if(me.eq.king.or..not.out1file)
write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
& (d_t(j,i+nres),j=1,3)
enddo
+ endif
+ if (lang.eq.0) then
c Zeroing the total angular momentum of the system
write(iout,*) "Calling the zero-angular
& momentum subroutine"
- endif
call inertia_tensor
c Getting the potential energy and forces and velocities and accelerations
call vcm_vel(vcm)
if(me.eq.king.or..not.out1file)then
write (iout,*) "vcm right after adjustment:"
write (iout,*) (vcm(j),j=1,3)
+ write (iout,*) "Initial velocities after adjustment"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
+ & (d_t(j,i+nres),j=1,3)
+ enddo
+ call flush(iout)
endif
- call flush(iout)
- write (iout,*) "init_MD before initial structure REST ",rest
- if (.not.rest) then
+ endif
+c write (iout,*) "init_MD before initial structure REST ",rest
+ if(start_from_model .and. (me.eq.king .or. .not. out1file))
+ & write(iout,*) 'START_FROM_MODELS is ON'
+ if(start_from_model .and. rest) then
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*)
+ & 'START_FROM_MODELS is OFF because the run is restarted'
+ write(iout,*) 'Remove restart keyword from input'
+ endif
+ endif
+c write (iout,*) "rest ",rest," start_from_model",start_from_model,
+c & " nmodel_start",nmodel_start," preminim",preminim
+ if (.not.rest) then
+ 122 continue
if (iranconf.ne.0) then
c 8/22/17 AL Loop to produce a low-energy random conformation
do iranmin=1,10
44 continue
else if (preminim) then
if (start_from_model) then
- i_model=iran_num(1,constr_homology)
- write (iout,*) 'starting from model ',i_model
- do i=1,2*nres
- do j=1,3
- c(j,i)=chomo(j,i,i_model)
+ n_model_try=0
+ fail=.true.
+ list_model_try=0
+ do while (fail .and. n_model_try.lt.nmodel_start)
+ write (iout,*) "n_model_try",n_model_try
+ do
+ i_model=iran_num(1,nmodel_start)
+ do k=1,n_model_try
+ if (i_model.eq.list_model_try(k)) exit
+ enddo
+ if (k.gt.n_model_try) exit
enddo
- enddo
- call int_from_cart(.true.,.false.)
- call sc_loc_geom(.false.)
- 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)
+ n_model_try=n_model_try+1
+ list_model_try(n_model_try)=i_model
+ if (me.eq.king .or. .not. out1file)
+ & write (iout,*) 'Trying to start from model ',
+ & pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=chomo(j,i,i_model)
+ enddo
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)
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ 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
- enddo
- endif
+ 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
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies before removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
! Remove SC overlaps if requested
- if (overlapsc) then
- write (iout,*) 'Calling OVERLAP_SC'
- call overlap_sc(fail)
- endif
+ if (overlapsc) then
+ write (iout,*) 'Calling OVERLAP_SC'
+ call overlap_sc(fail)
+ if (fail) then
+ write (iout,*)
+ & "Failed to remove overlap from model",i_model
+ cycle
+ endif
+ endif
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies after removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
+#ifdef SEARCHSC
! Search for better SC rotamers if requested
- if (searchsc) then
- call sc_move(2,nres-1,10,1d10,nft_sc,etot)
- print *,'SC_move',nft_sc,etot
- if (me.eq.king.or..not.out1file)
- & write(iout,*) 'SC_move',nft_sc,etot
+ if (searchsc) then
+ call sc_move(2,nres-1,10,1d10,nft_sc,etot)
+ print *,'SC_move',nft_sc,etot
+ if (me.eq.king.or..not.out1file)
+ & write(iout,*) 'SC_move',nft_sc,etot
+ endif
+ call etotal(energia(0))
+#endif
+ enddo
+ call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),
+ & 1,MPI_INTEGER,king,CG_COMM,IERROR)
+ if (n_model_try.gt.nmodel_start .and.
+ & (me.eq.king .or. out1file)) then
+ write (iout,*)
+ & "All models have irreparable overlaps. Trying randoms starts."
+ iranconf=1
+ i_model=nmodel_start+1
+ goto 122
+ endif
+ else
+! Remove SC overlaps if requested
+ if (overlapsc) then
+ write (iout,*) 'Calling OVERLAP_SC'
+ call overlap_sc(fail)
+ if (fail) then
+ write (iout,*)
+ & "Failed to remove overlap"
+ endif
+ endif
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies after removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
endif
- call etotal(energia(0))
C 8/22/17 AL Minimize initial structure
if (dccart) then
if (me.eq.king.or..not.out1file) write(iout,*)
- & 'Minimizing initial PDB structure: Calling MINIM_DC'
+ & 'Minimizing initial PDB structure: Calling MINIM_DC'
call minim_dc(etot,iretcode,nfun)
else
call geom_to_var(nvar,varia)
if(me.eq.king.or..not.out1file) write (iout,*)
- & 'Minimizing initial PDB structure: Calling MINIMIZE.'
+ & 'Minimizing initial PDB structure: Calling MINIMIZE.'
call minimize(etot,varia,iretcode,nfun)
call var_to_geom(nvar,varia)
- endif
- if (me.eq.king.or..not.out1file)
+#ifdef LBFGS
+ if (me.eq.king.or..not.out1file)
+ & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+#else
+ if (me.eq.king.or..not.out1file)
& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
- if(me.eq.king.or..not.out1file)
+ if(me.eq.king.or..not.out1file)
& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+#endif
+ endif
endif
- endif
+ if (nmodel_start.gt.0 .and. me.eq.king) then
+ write (iout,'(a)') "Task Starting model"
+ do i=0,nodes-1
+ if (i_start_models(i).gt.nmodel_start) then
+ write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
+ else
+ write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i))
+ & (:ilen(pdbfiles_chomo(i_start_models(i))))
+ endif
+ enddo
+ endif
+ endif ! .not. rest
call chainbuild_cart
call kinetic(EK)
if (tbf) then
call verlet_bath
- endif
+ endif
kinetic_T=2.0d0/(dimen3*Rb)*EK
+ write (iout,*) "Initial kinetic energy",EK," kinetic T",kinetic_T
if(me.eq.king.or..not.out1file)then
call cartprint
call intout
c write (iout,*) "PotE-homology",potE
totE=EK+potE
itime=0
+ itime_mat=itime
if (ntwe.ne.0) call statout(itime)
if(me.eq.king.or..not.out1file)
& write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
integer i,ii,j,k,l,ind
double precision anorm_distr
- logical lprn /.true./
+ logical lprn /.false./
#ifdef FIVEDIAG
integer ichain,n,innt,inct,ibeg,ierr
double precision work(8*maxres6)
!
#define WLOS
#ifdef WLOS
+ if (nnt.eq.1) then
+ d_t(:,0)=d_t(:,1)
+ endif
do i=1,nres
if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
do j=1,3
enddo
end if
enddo
+ d_t(:,nres)=0.0d0
d_t(:,nct)=0.0d0
+ d_t(:,2*nres)=0.0d0
+ if (nnt.gt.1) then
+ d_t(:,0)=d_t(:,1)
+ d_t(:,1)=0.0d0
+ endif
c d_a(:,0)=d_a(:,1)
c d_a(:,1)=0.0d0
c write (iout,*) "Shifting accelerations"
- do ichain=1,nchain
+ do ichain=2,nchain
c write (iout,*) "ichain",chain_border1(1,ichain)-1,
c & chain_border1(1,ichain)
d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
& d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
d_t(:,chain_border(2,ichain-1))=0.0d0
enddo
+ do ichain=2,nchain
+ write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+ & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=
+ & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
#else
ibeg=0
c do j=1,3
dc(j,0)=dc_work(j)
d_t(j,0)=d_t_work(j)
enddo
- ind=3
- do i=nnt,nct-1
+ ind=3
+ do i=nnt,nct-1
do j=1,3
dc(j,i)=dc_work(ind+j)
d_t(j,i)=d_t_work(ind+j)
projects / unres.git / blobdiff