Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / unres / MD.F90
index 2b412fd..52d7ec3 100644 (file)
       character(len=50) :: tytul
 !el      common /gucio/ cm
       integer :: i,j,nharp
-      integer,dimension(4,nres/3) :: iharp     !(4,nres/3)(4,maxres/3)
+      integer,dimension(4,nres) :: iharp       !(4,nres/3)(4,maxres/3)
 !      logical :: ovrtim
       real(kind=8) :: tt0,scalfac
       integer :: nres2,itime
           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)
        write (iout,*) "vcm right after adjustment:"
        write (iout,*) (vcm(j),j=1,3) 
       endif
+
+
+  
+!         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'
             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
 !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