Adam's unres update
[unres.git] / source / unres / src-HCD-5D / unres.F
index 76bd280..2e0ebaf 100644 (file)
@@ -20,7 +20,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       include 'COMMON.GEO'
       include 'COMMON.HEADER'
       include 'COMMON.CONTROL'
-      include 'COMMON.CONTACTS'
+c      include 'COMMON.CONTACTS'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
@@ -59,7 +59,8 @@ c      call memmon_print_usage()
 C Read force field parameters and job setup data
       call readrtns
 C
-c      write (iout,*) "After readrtns"
+      write (iout,*) "After readrtns"
+      call flush(iout)
       call cartprint
       call intout
       if (me.eq.king .or. .not. out1file) then
@@ -73,8 +74,8 @@ c      write (iout,*) "After readrtns"
       call flush(iout)
 C
       if (modecalc.eq.-2) then
-        call test
-        stop
+c        call test
+c        stop
       else if (modecalc.eq.-1) then
         write(iout,*) "call check_sc_map next"
         call check_bond
@@ -86,9 +87,14 @@ C Fine-grain slaves just do energy and gradient components.
         call ergastulum ! slave workhouse in Latin
       else
 #endif
+      if (indpdb.eq.0 .and. .not.read_cart) call chainbuild
+      if (indpdb.ne.0 .or. read_cart) then
+        dc(1,0)=c(1,1)
+        dc(2,0)=c(2,1)
+        dc(3,0)=c(3,1)
+      endif
       if (modecalc.eq.0) then
         write (iout,*) "Calling exec_eeval_or_minim"
-        call cartprint
         call exec_eeval_or_minim
       else if (modecalc.eq.1) then
         call exec_regularize
@@ -199,7 +205,7 @@ c---------------------------------------------------------------------------
       include 'COMMON.GEO'
       include 'COMMON.HEADER'
       include 'COMMON.CONTROL'
-      include 'COMMON.CONTACTS'
+c      include 'COMMON.CONTACTS'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
@@ -207,19 +213,24 @@ c---------------------------------------------------------------------------
       include 'COMMON.REMD'
       include 'COMMON.MD'
       include 'COMMON.SBRIDGE'
-      integer i,icall,iretcode,nfun
+      integer i,it,icall,iretcode,nfun
       common /srutu/ icall
       integer nharp,iharp(4,maxres/3)
       integer nft_sc
-      logical fail
+      logical fail,secondary_str /.true./
       double precision energy(0:n_ene),etot,etota
       double precision energy_long(0:n_ene),energy_short(0:n_ene)
       double precision rms,frac,frac_nn,co
       double precision varia(maxvar)
       double precision time00,time1,time_ene,evals
+#ifdef LBFGS
+      character*9 status
+      integer niter
+      common /lbfgstat/ status,niter,nfun
+#endif
       integer ilen
-      if (indpdb.eq.0)     call chainbuild
-      if (indpdb.ne.0) then
+      if (indpdb.eq.0 .and. .not.read_cart) call chainbuild
+      if (indpdb.ne.0 .or. read_cart) then
       dc(1,0)=c(1,1)
       dc(2,0)=c(2,1)
       dc(3,0)=c(3,1)
@@ -232,6 +243,24 @@ c---------------------------------------------------------------------------
       write (iout,*) "Energy evaluation/minimization"
       call chainbuild_cart
 c      print *,'dc',dc(1,0),dc(2,0),dc(3,0)
+      if (nran_start.gt.0) then
+        write (iout,*) 
+     & "Chains will be regenerated starting from residue",nran_start
+        do it=1,100
+        call gen_rand_conf_mchain(nran_start,*10)
+        write (iout,*) "Conformation successfully generated",it
+        goto 11
+   10   write (iout,*) "Problems with regenerating chains",it
+        enddo
+   11   continue
+        write (iout,*) "Cartesian coords after chain rebuild"
+        call cartprint
+        call chainbuild_cart
+        write (iout,*) "Cartesian coords after chainbuild_ecart"
+        call cartprint
+        call int_from_cart1(.false.)
+        call intout
+      endif
       if (split_ene) then
        print *,"Processor",myrank," after chainbuild"
        icall=1
@@ -263,16 +292,29 @@ c      print *,"after etotal"
       etota = energy(0)
       etot =etota
       call enerprint(energy(0))
+      if (secondary_str) then
       call hairpin(.true.,nharp,iharp)
 c        print *,'after hairpin'
       call secondary2(.true.)
+      endif
 c        print *,'after secondary'
       if (minim) then
 crc overlap test
+        if (indpdb.ne.0 .and. .not.dccart) then 
+          call bond_regular
+          call chainbuild_extconf
+          call etotal(energy(0))
+          write (iout,*) "After bond regularization"
+          call enerprint(energy(0))
+        endif
+
         if (overlapsc) then 
-          print *, 'Calling OVERLAP_SC'
+          write (iout,*) 'Calling OVERLAP_SC'
           call overlap_sc(fail)
-          print *,"After overlap_sc"
+          write (iout,*) "After overlap_sc"
+c        cartname=prefix(:ilen(prefix))//'.x'
+c        potE=etot
+c        call cartoutx(0.0d0)
         endif 
 
         if (searchsc) then 
@@ -290,12 +332,8 @@ crc overlap test
 #endif
           call minim_dc(etot,iretcode,nfun)
         else
-          if (indpdb.ne.0) then 
-            call bond_regular
-            call chainbuild_extconf
-          endif
           call geom_to_var(nvar,varia)
-          print *,'Calling MINIMIZE.'
+c          print *,'Calling MINIMIZE.'
 #ifdef MPI
           time1=MPI_WTIME()
 #else
@@ -303,7 +341,11 @@ crc overlap test
 #endif
           call minimize(etot,varia,iretcode,nfun)
         endif
+#ifdef LBFGS
+        print *,'LBFGS return code is',status,' eval ',nfun
+#else
         print *,'SUMSL return code is',iretcode,' eval ',nfun
+#endif
 #ifdef MPI
         evals=nfun/(MPI_WTIME()-time1)
 #else
@@ -311,29 +353,32 @@ crc overlap test
 #endif
         print *,'# eval/s',evals
         print *,'refstr=',refstr
+        if (secondary_str) then
         call hairpin(.false.,nharp,iharp)
-        print *,'after hairpin'
+c        print *,'after hairpin'
         call secondary2(.true.)
-        print *,'after secondary'
+c        print *,'after secondary'
+        endif
         call etotal(energy(0))
         etot = energy(0)
         call enerprint(energy(0))
-
         call intout
-        if (out_int) call briefout(0,etot)
-        if (out_cart) then
-          cartname=prefix(:ilen(prefix))//'.x'
-          potE=etot
-          call cartoutx(0.0d0)
-        endif
-        if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+#ifdef LBFGS
+          write (iout,'(a,a9)') 'LBFGS return code:',status
+          write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
+          write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
+#else
           write (iout,'(a,i3)') 'SUMSL return code:',iretcode
           write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
           write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
-      else
-        print *,'refstr=',refstr
-        if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
-        call briefout(0,etot)
+#endif
+      endif
+      if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+      if (out_int) call briefout(0,etot)
+      if (out_cart) then
+        cartname=prefix(:ilen(prefix))//'.x'
+        potE=etot
+        call cartoutx(0.0d0)
       endif
       if (outpdb) call pdbout(etot,titel(:50),ipdb)
       if (outmol2) call mol2out(etot,titel(:32))
@@ -353,7 +398,7 @@ c---------------------------------------------------------------------------
       include 'COMMON.GEO'
       include 'COMMON.HEADER'
       include 'COMMON.CONTROL'
-      include 'COMMON.CONTACTS'
+c      include 'COMMON.CONTACTS'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
@@ -364,6 +409,11 @@ c---------------------------------------------------------------------------
       double precision energy(0:n_ene)
       double precision etot,rms,frac,frac_nn,co
       integer iretcode
+#ifdef LBFGS
+      character*9 status
+      integer niter,nfun
+      common /lbfgstat/ status,niter,nfun
+#endif
 
       call gen_dist_constr
       call sc_conf
@@ -378,7 +428,11 @@ c---------------------------------------------------------------------------
       if (outpdb) call pdbout(etot,titel(:50),ipdb)
       if (outmol2) call mol2out(etot,titel(:32))
       if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+#ifdef LBFGS
+      write (iout,'(a,a9)') 'LBFGS return code:',status
+#else
       write (iout,'(a,i3)') 'SUMSL return code:',iretcode
+#endif
       return
       end
 c---------------------------------------------------------------------------
@@ -440,7 +494,7 @@ c---------------------------------------------------------------------------
       include 'COMMON.GEO'
       include 'COMMON.HEADER'
       include 'COMMON.CONTROL'
-      include 'COMMON.CONTACTS'
+c      include 'COMMON.CONTACTS'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
@@ -739,7 +793,7 @@ c---------------------------------------------------------------------------
       include 'COMMON.GEO'
       include 'COMMON.HEADER'
       include 'COMMON.CONTROL'
-      include 'COMMON.CONTACTS'
+c      include 'COMMON.CONTACTS'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
@@ -825,7 +879,8 @@ c      enddo
       goto (10,20,30) icheckgrad
   10  call check_ecartint
       return
-  20  call check_cartgrad
+  20  write (iout,*) 
+     & "Checking the gradient of Cartesian coordinates disabled."
       return
   30  call check_eint
       return