rescore=0.2 rescore2=1.0 mdpdb append
[unres.git] / source / unres / src_MD-M / MD_A-MTS.F
index 7a9e307..d75964a 100644 (file)
@@ -36,6 +36,7 @@ c------------------------------------------------
       common /gucio/ cm
       integer itime
       logical ovrtim
+      integer nharp,iharp(4,maxres/3)
 c
 #ifdef MPI
       if (ilen(tmpdir).gt.0)
@@ -51,6 +52,33 @@ c
       t_enegrad=0.0d0
       t_sdsetup=0.0d0
       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
+      write (iout,'(/a)')
+     &  "Cartesian coordinates of the initial 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
+      write (iout,'(/a)') 
+     &  "Initial dC vectors of the chain"
+      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,(dc(j,ires),j=1,3),
+     &    (dc(j,ires+nres),j=1,3)
+      enddo
+      write (iout,'(/a)') 
+     &  "Initial dC_norm vectors of the chain"
+      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,(dc_norm(j,ires),j=1,3),
+     &    (dc_norm(j,ires+nres),j=1,3)
+      enddo
 #ifdef MPI
       tt0=MPI_Wtime()
 #else
@@ -1443,7 +1471,7 @@ c  Set up the initial conditions of a MD simulation
       include 'COMMON.NAMES'
       include 'COMMON.REMD'
       real*8 energia_long(0:n_ene),
-     &  energia_short(0:n_ene),vcm(3),incr(3)
+     &  energia_short(0:n_ene),energia(0:n_ene),vcm(3),incr(3)
       double precision cm(3),L(3),xv,sigv,lowb,highb
       double precision varia(maxvar)
       character*256 qstr
@@ -1471,9 +1499,15 @@ c Open the pdb file for snapshotshots
         if (ilen(tmpdir).gt.0) 
      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
      &      liczba(:ilen(liczba))//".pdb")
+#if defined(AIX) || defined(PGI)
         open(ipdb,
      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
-     &  //".pdb")
+     &  //".pdb", position='append')
+#else
+        open(ipdb,
+     &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
+     &  //".pdb", access='append')
+#endif
       else
 #ifdef NOXDR
         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) 
@@ -1609,11 +1643,38 @@ c Removing the velocity of the center of mass
        write (iout,*) (vcm(j),j=1,3) 
       endif
       if (.not.rest) then              
-         call chainbuild
+         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)
+           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)
+            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
+          call etotal(energia(0))
+         endif
+c         call chainbuild_cart
+         write (iout,*) "Initial energies"
+         call enerprint(energia(0))
          write (iout,*) "PREMINIM ",preminim
          if(iranconf.ne.0 .or. preminim) then
           if (overlapsc) then 
-           print *, 'Calling OVERLAP_SC'
+           write (iout,*) 'Calling OVERLAP_SC'
            call overlap_sc(fail)
           endif 
 
@@ -1640,6 +1701,81 @@ c Removing the velocity of the center of mass
              call enerprint(potEcomp)
           endif
          endif
+         if(iranconf.ne.0) then
+c 8/22/17 AL Loop to produce a low-energy random conformation
+          DO iranmin=1,10
+          if (overlapsc) then
+           if(me.eq.king.or..not.out1file) 
+     &        write (iout,*) 'Calling OVERLAP_SC'
+           call overlap_sc(fail)
+          endif
+
+          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
+
+          if(dccart)then
+           print *, 'Calling MINIM_DC'
+           call minim_dc(etot,iretcode,nfun)
+          else
+           call geom_to_var(nvar,varia)
+           print *,'Calling MINIMIZE.'
+           call minimize(etot,varia,iretcode,nfun)
+           call var_to_geom(nvar,varia)
+          endif
+          if(me.eq.king.or..not.out1file)
+     &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+
+          if (isnan(etot) .or. etot.gt.1.0d4) then
+            write (iout,*) "Energy too large",etot,
+     &        " trying another random conformation"
+            do itrial=1,100
+              itmp=1
+              call gen_rand_conf(itmp,*30)
+              goto 40
+   30         write (iout,*) 'Failed to generate random conformation',
+     &          ', itrial=',itrial
+              write (*,*) 'Processor:',me,
+     &          ' Failed to generate random conformation',
+     &          ' itrial=',itrial
+              call intout
+#ifdef AIX
+              call flush_(iout)
+#else
+              call flush(iout)
+#endif
+            enddo
+            write (iout,'(a,i3,a)') 'Processor:',me,
+     &        ' error in generating random conformation.'
+            write (*,'(a,i3,a)') 'Processor:',me,
+     &        ' error in generating random conformation.'
+            call flush(iout)
+#ifdef MPI
+            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+#else
+            stop
+#endif
+   40       continue
+          else
+            goto 44
+          endif
+          ENDDO
+
+          write (iout,'(a,i3,a)') 'Processor:',me,
+     &        ' failed to generate a low-energy random conformation.'
+            write (*,'(a,i3,a)') 'Processor:',me,
+     &        ' failed to generate a low-energy random conformation.'
+            call flush(iout)
+#ifdef MPI
+            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+#else
+            stop
+#endif
+   44     continue
+         endif
       endif      
       call chainbuild_cart
       call kinetic(EK)
@@ -1676,11 +1812,11 @@ c Removing the velocity of the center of mass
      &   "Time step reduced to",d_time,
      &   " because of too large initial acceleration."
       endif
-C      if(me.eq.king.or..not.out1file)then 
-C       write(iout,*) "Potential energy and its components"
-C       call enerprint(potEcomp)
+      if(me.eq.king.or..not.out1file)then 
+       write(iout,*) "Potential energy and its components"
+       call enerprint(potEcomp)
 c       write(iout,*) (potEcomp(i),i=0,n_ene)
-C      endif
+      endif
       potE=potEcomp(0)-potEcomp(20)
       totE=EK+potE
       itime=0