remd with homol_nset>1 working with FG parallelization
[unres.git] / source / unres / src_MD / readrtns.F
index 056bfd4..69c7fe8 100644 (file)
@@ -8,6 +8,7 @@
       include 'COMMON.CONTROL'
       include 'COMMON.SBRIDGE'
       include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
       logical file_exist
 C Read force-field parameters except weights
       call parmread
@@ -130,6 +131,7 @@ C Set up the time limit (caution! The time must be input in minutes!)
       sideadd=(index(controlcard,'SIDEADD').gt.0)
       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
       outpdb=(index(controlcard,'PDBOUT').gt.0)
+      outx=(index(controlcard,'XOUT').gt.0)
       outmol2=(index(controlcard,'MOL2OUT').gt.0)
       pdbref=(index(controlcard,'PDBREF').gt.0)
       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
@@ -386,6 +388,11 @@ C
       large = index(controlcard,"LARGE").gt.0
       print_compon = index(controlcard,"PRINT_COMPON").gt.0
       rattle = index(controlcard,"RATTLE").gt.0
+      preminim = index(controlcard,"PREMINIM").gt.0
+      if (preminim) then
+        dccart=(index(controlcard,'CART').gt.0)
+        call read_minim
+      endif
 c  if performing umbrella sampling, fragments constrained are read from the fragment file 
       nset=0
       if(usampl) then
@@ -426,6 +433,10 @@ c  if performing umbrella sampling, fragments constrained are read from the frag
        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
        if (rattle) write (iout,'(a60)') 
      &  "Rattle algorithm used to constrain the virtual bonds"
+       if (preminim .or. iranconf.gt.0) then
+         write (iout,'(a60)')
+     &      "Initial structure will be energy-minimized" 
+       endif
       endif
       reset_fricmat=1000
       if (lang.gt.0) then
@@ -627,6 +638,7 @@ C
       common /pizda/ itype_pdb
       logical seq_comp,fail
       double precision energia(0:n_ene)
+      
       integer ilen
       external ilen
 C
@@ -922,6 +934,17 @@ C 12/1/95 Added weight for the multi-body term WCORR
   34    continue
 c        print *,'Begin reading pdb data'
         call readpdb
+        do i=1,2*nres
+          do j=1,3
+            crefjlee(j,i)=c(j,i)
+          enddo
+        enddo
+#ifdef DEBUG
+        do i=1,nres
+          write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
+     &      (crefjlee(j,i+nres),j=1,3)
+        enddo
+#endif
 c        print *,'Finished reading pdb data'
         if(me.eq.king.or..not.out1file)
      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
@@ -1165,6 +1188,44 @@ c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 
       if (constr_homology.gt.0) then
         call read_constr_homology
+        if (indpdb.gt.0 .or. pdbref) then
+          do i=1,2*nres
+            do j=1,3
+              c(j,i)=crefjlee(j,i)
+              cref(j,i)=crefjlee(j,i)
+            enddo
+          enddo
+        endif
+#ifdef DEBUG
+        write (iout,*) "Array C"
+        do i=1,nres
+          write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
+     &      (c(j,i+nres),j=1,3)
+        enddo
+        write (iout,*) "Array Cref"
+        do i=1,nres
+          write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
+     &      (cref(j,i+nres),j=1,3)
+        enddo
+#endif
+       call int_from_cart1(.false.)
+       call sc_loc_geom(.false.)
+       do i=1,nres
+         thetaref(i)=theta(i)
+         phiref(i)=phi(i)
+       enddo
+       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
       else
         homol_nset=0
       endif
@@ -1182,6 +1243,8 @@ C initial geometry.
           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
      &     write (iout,'(a)') 'Initial geometry will be read in.'
           if (read_cart) then
+            read (inp,*) time,potE,uconst,t_bath,
+     &       nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
             read(inp,'(8f10.5)',end=36,err=36)
      &       ((c(l,k),l=1,3),k=1,nres),
      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
@@ -1200,7 +1263,7 @@ C initial geometry.
                 enddo
               endif
             enddo
-            return
+c            return
           else
             call read_angles(inp,*36)
           endif
@@ -1284,6 +1347,8 @@ C Generate distance constraints, if the PDB structure is to be regularized.
       if (nthread.gt.0) then
         call read_threadbase
       endif
+      write (iout,*) "READRTNS: Calling setup_var"
+      call flush(iout)
       call setup_var
       if (me.eq.king .or. .not. out1file)
      & call intout
@@ -2262,6 +2327,8 @@ c      print *,"Processor",myrank," fg_rank",fg_rank
      &  //'.pdb'
       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
      &  liczba(:ilen(liczba))//'.mol2'
+      cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
+     &  liczba(:ilen(liczba))//'.x'
       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
      &  liczba(:ilen(liczba))//'.stat'
       if (lentmp.gt.0)
@@ -2280,6 +2347,7 @@ c      print *,"Processor",myrank," fg_rank",fg_rank
       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
+      cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
       if (lentmp.gt.0)
      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
@@ -2440,33 +2508,36 @@ c-------------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.MD'
       include 'COMMON.CONTROL'
+      integer iset1
       read(inp,*) nset,nfrag,npair,nfrag_back
       if(me.eq.king.or..not.out1file)
      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
      &  " nfrag_back",nfrag_back
-      do iset=1,nset
-         read(inp,*) mset(iset)
+      do iset1=1,nset
+         read(inp,*) mset(iset1)
        do i=1,nfrag
-         read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset), 
-     &     qinfrag(i,iset)
+         read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1), 
+     &     qinfrag(i,iset1)
          if(me.eq.king.or..not.out1file)
-     &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
-     &     ifrag(2,i,iset), qinfrag(i,iset)
+     &    write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
+     &     ifrag(2,i,iset1), qinfrag(i,iset1)
        enddo
        do i=1,npair
-        read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset), 
-     &    qinpair(i,iset)
+        read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1), 
+     &    qinpair(i,iset1)
         if(me.eq.king.or..not.out1file)
-     &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
-     &    ipair(2,i,iset), qinpair(i,iset)
+     &   write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
+     &    ipair(2,i,iset1), qinpair(i,iset1)
        enddo 
        do i=1,nfrag_back
-        read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
-     &     wfrag_back(3,i,iset),
-     &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+        read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
+     &     wfrag_back(3,i,iset1),
+     &     ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
         if(me.eq.king.or..not.out1file)
-     &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
-     &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+     &   write(iout,*) "A",i,wfrag_back(1,i,iset1),
+     &   wfrag_back(2,i,iset1),
+     &   wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
+     &   ifrag_back(2,i,iset1)
        enddo 
       enddo
       return
@@ -2672,22 +2743,21 @@ c Alternative: reading from input
       call readi(controlcard,"HOMOL_NSET",homol_nset,1)       
       if (homol_nset.gt.1)then
          call card_concat(controlcard)
-         read(controlcard,*) (waga_dist1(i),i=1,homol_nset) 
-         call card_concat(controlcard)
-         read(controlcard,*) (waga_angle1(i),i=1,homol_nset) 
-         write(iout,*) "iset distance_weight angle_weight"
-         do i=1,homol_nset
-           write(iout,*) i,waga_dist1(i),waga_angle1(i)
-         enddo
+         read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
+         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+          write(iout,*) "iset homology_weight "
+          do i=1,homol_nset
+           write(iout,*) i,waga_homology(i)
+          enddo
+         endif
          iset=mod(kolor,homol_nset)+1
       else
        iset=1
-c       call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
-c       call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
+       waga_homology(1)=1.0
       endif
 
-      write (iout,*) "nnt",nnt," nct",nct
-      call flush(iout)
+cd      write (iout,*) "nnt",nnt," nct",nct
+cd      call flush(iout)
 
 
       lim_odl=0
@@ -2734,6 +2804,7 @@ c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
 
+        unres_pdb=.false.
         call readpdb
 c
 c     Distance restraints
@@ -2773,7 +2844,7 @@ c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
 c         enddo
 c 1401   continue
 c         close (ientin)
-        if (waga_dist.gt.0.0d0) then
+        if (waga_dist.ne.0.0d0) then
           ii=0
           do i = nnt,nct-2 ! right? without parallel.
             do j=i+2,nct ! right?
@@ -2938,40 +3009,41 @@ c              read (ientin,*) sigma_d(k,i) ! 1st variant
         endif
         close(ientin)
       enddo
-      if (waga_dist.gt.0.0d0) lim_odl=ii
+      if (waga_dist.ne.0.0d0) lim_odl=ii
       if (constr_homology.gt.0) call homology_partition
       if (constr_homology.gt.0) call init_int_table
-      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
-     & "lim_xx=",lim_xx
+cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
+cd     & "lim_xx=",lim_xx
 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 c
 c Print restraints
 c
       if (.not.lprn) return
-      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
-      write (iout,*) "Distance restraints from templates"
-      do ii=1,lim_odl
+cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+       write (iout,*) "Distance restraints from templates"
+       do ii=1,lim_odl
        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
-      enddo
-      write (iout,*) "Dihedral angle restraints from templates"
-      do i=nnt+3,lim_dih
+       enddo
+       write (iout,*) "Dihedral angle restraints from templates"
+       do i=nnt+3,lim_dih
         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
-      enddo
-      write (iout,*) "Virtual-bond angle restraints from templates"
-      do i=nnt+2,lim_theta
+       enddo
+       write (iout,*) "Virtual-bond angle restraints from templates"
+       do i=nnt+2,lim_theta
         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
-      enddo
-      write (iout,*) "SC restraints from templates"
-      do i=nnt,lim_xx
+       enddo
+       write (iout,*) "SC restraints from templates"
+       do i=nnt,lim_xx
         write(iout,'(i5,10(4f8.2,4x))') i,
      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
-      enddo
-
+       enddo
+      endif
 c -----------------------------------------------------------------
       return
       end