!el external ilen
character(len=50) :: tytul
!el common /gucio/ cm
- integer :: itime,i,j,nharp
- integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
+ integer :: i,j,nharp
+ integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
! logical :: ovrtim
real(kind=8) :: tt0,scalfac
- integer :: nres2
+ integer :: nres2,itime
nres2=2*nres
print *, "ENTER MD"
!
stop
#endif
endif
+ itime_mat=itime
if (ntwe.ne.0) then
if (mod(itime,ntwe).eq.0) then
call statout(itime)
use comm_gucio
use control, only:tcpu
use control_data
+ use minimm, only:minim_dc
#ifdef MPI
include 'mpif.h'
integer :: ierror,ierrcode
integer :: count,rstcount !ilen,
!el external ilen
character(len=50) :: tytul
- integer :: maxcount_scale = 20
+ integer :: maxcount_scale = 30
!el common /gucio/ cm
!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
!el common /stochcalc/ stochforcvec
- integer :: itime,icount_scale,itime_scal,i,j,ifac_time
- logical :: scale
+ integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
+ logical :: scalel
real(kind=8) :: epdrift,tt0,fac_time
!
if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
- scale=.true.
+ scalel=.true.
icount_scale=0
if (lang.eq.1) then
call sddir_precalc
#endif
endif
itime_scal=0
- do while (scale)
+ do while (scalel)
icount_scale=icount_scale+1
+! write(iout,*) "icount_scale",icount_scale,scalel
if (icount_scale.gt.maxcount_scale) then
write (iout,*) &
"ERROR: too many attempts at scaling down the time step. ",&
call etotal(potEcomp)
! AL 4/17/17: Reduce the steps if NaNs occurred.
if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
- d_time=d_time/2
+ call enerprint(potEcomp)
+ d_time=d_time/10.0
+ if (icount_scale.gt.15) then
+ write (iout,*) "Tu jest problem",potEcomp(0),d_time
+! call gen_rand_conf(1,*335)
+! call minim_dc(potEcomp(0),iretcode,100)
+
+! call zerograd
+! call etotal(potEcomp)
+! write(iout,*) "needed to repara,",potEcomp
+ endif
cycle
+! 335 write(iout,*) "Failed genrand"
+! cycle
endif
! end change
if (large.and. mod(itime,ntwe).eq.0) &
call max_accel
amax=amax/(itime_scal+1)**2
call predict_edrift(epdrift)
+! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
+ scalel=.false.
+! write (iout,*) "before enter if",scalel,icount_scale
if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
+! write(iout,*) "I enter if"
! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
- scale=.true.
+ scalel=.true.
ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
/dlog(2.0d0)+1
itime_scal=itime_scal+ifac_time
endif
#endif
endif
- scale=.false.
endif
enddo
! Calculate the kinetic and the total energy and the kinetic temperature
!el common /gucio/ cm
!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
!el common /stochcalc/ stochforcvec
- integer :: itime,itt,i,j,itsplit
+ integer :: itt,i,j,itsplit,itime
logical :: scale
!el common /cipiszcze/ itt
do j=1,3
adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
+! write(iout,*) "adt",adt,"ads",adt2
dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
d_t(j,inres)=d_t_old(j,inres)+adt
use minimm, only:minim_dc,minimize,sc_move
use io_config, only:readrst
use io, only:statout
+ use random, only: iran_num
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
#ifdef MP
character(len=16) :: form
integer :: IERROR,ERRCODE
#endif
- integer :: iranmin,itrial,itmp
+ integer :: iranmin,itrial,itmp,n_model_try,k, &
+ i_model
+ integer, dimension(:),allocatable :: list_model_try
+ integer, dimension(0:nodes-1) :: i_start_models
! include 'COMMON.SETUP'
! include 'COMMON.CONTROL'
! include 'COMMON.VAR'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
! include 'COMMON.REMD'
- real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
+ real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
real(kind=8),dimension(3) :: vcm,incr,L
real(kind=8) :: xv,sigv,lowb,highb
real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
character(len=50) :: tytul
logical :: file_exist
!el common /gucio/ cm
- integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
+ integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
real(kind=8) :: etot,tt0
logical :: fail
*dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
enddo
endif
+
! Open the pdb file for snapshotshots
#ifdef MPI
if(mdpdb) then
write (iout,*) "vcm right after adjustment:"
write (iout,*) (vcm(j),j=1,3)
endif
- if (.not.rest) then
+
+
+
! call chainbuild
+
+ if ((.not.rest).or.(forceminim)) then
+ if (forceminim) call chainbuild_cart
+ 122 continue
if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
if (overlapsc) then
print *, 'Calling OVERLAP_SC'
if(dccart)then
print *, 'Calling MINIM_DC'
call minim_dc(etot,iretcode,nfun)
+ call int_from_cart1(.false.)
else
call geom_to_var(nvar,varia)
print *,'Calling MINIMIZE.'
stop
#endif
44 continue
+ else if (preminim) then
+ if (start_from_model) then
+ 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
+ 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
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ dc(:,0)=c(:,1)
+ 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
+ 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)
+ 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
+ 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
+! 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'
+ 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.'
+ call minimize(etot,varia,iretcode,nfun)
+ call var_to_geom(nvar,varia)
+#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)&
+ write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+#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
call chainbuild_cart
potE=potEcomp(0)-potEcomp(20)
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:", &
!el ginvfric(2*nres,2*nres) !maxres2=2*maxres
!el common /przechowalnia/ ginvfric
- logical :: lprn = .false., checkmode = .false.
+ logical :: lprn, checkmode
integer :: i,j,ind,k,nres2,nres6,mnum
nres2=2*nres
nres6=6*nres
-
+ lprn=.false.
+ checkmode=.false.
+! if (large) lprn=.true.
+! if (large) checkmode=.true.
if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
do i=0,nres2
projects / unres4.git / blobdiff