updates from src_MD to src_MD-M for unres

diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN
index 84ebf1b..dd53a8e 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CHAIN
+++ b/source/unres/src_MD-M/COMMON.CHAIN
@@ -3,7 +3,7 @@
      &  afmbeg
       double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r,
      & prod,rt,dc_work,cref,crefjlee,chain_rep,dc_norm2,velAFMconst,
-     & totTafm
+     & totTafm,chomo
       common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
      & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
      & dc_norm2(3,0:maxres2),
@@ -24,4 +24,4 @@
       common /afm/  forceAFMconst, distafminit,afmend,afmbeg,
      & velAFMconst,
      & totTafm
-
+      common /chomo_models/ chomo(3,maxres2+2,max_template)

diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL
index 0a868ce..35fb03d 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CONTROL
+++ b/source/unres/src_MD-M/COMMON.CONTROL
@@ -2,19 +2,21 @@
      & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
      & constr_homology,homol_nset
       real*8 waga_homology
-      real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut
+      real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut,
+     &  dist2_cut
       logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
      &                 sideadd,lsecondary,read_cart,unres_pdb,
      &                 vdisulf,searchsc,lmuca,dccart,extconf,out1file,
-     &                 gnorm_check,gradout,split_ene,read2sigma
+     &                 gnorm_check,gradout,split_ene,read2sigma,
+     &                 start_from_model
       common /cntrl/ modecalc,iscode,indpdb,indback,indphi,iranconf,
      & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint,
      & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
      & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
      & constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog,
      & selfguide,
-     & constr_homology,homol_nset,read2sigma
+     & constr_homology,homol_nset,read2sigma,start_from_model
       common /homol/ waga_homology(maxprocs/20),
-     & waga_dist, waga_angle, waga_theta, waga_d, dist_cut
+     & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut
 C... minim = .true. means DO minimization.
 C... energy_dec = .true. means print energy decomposition matrix

diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS
index 22c4844..70ce031 100644 (file)
--- a/source/unres/src_MD-M/DIMENSIONS
+++ b/source/unres/src_MD-M/DIMENSIONS
@@ -16,7 +16,7 @@ C Max. number of coarse-grain processors
       parameter (max_cg_procs=maxprocs)
 C Max. number of AA residues
       integer maxres
-      parameter (maxres=600)
+      parameter (maxres=800)
 C Appr. max. number of interaction sites
       integer maxres2,maxres6,mmaxres2
       parameter (maxres2=2*maxres,maxres6=6*maxres)

diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F
index 7a9e307..deef97c 100644 (file)
--- a/source/unres/src_MD-M/MD_A-MTS.F
+++ b/source/unres/src_MD-M/MD_A-MTS.F
@@ -1609,6 +1609,30 @@ c Removing the velocity of the center of mass
        write (iout,*) (vcm(j),j=1,3) 
       endif
       if (.not.rest) then              
+         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
+         endif
          call chainbuild
          write (iout,*) "PREMINIM ",preminim
          if(iranconf.ne.0 .or. preminim) then

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 79cd1b9..b3c4c22 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -6896,7 +6896,18 @@ c
            endif
          enddo
          
-         min_odl=minval(distancek)
+c         min_odl=minval(distancek)
+         do kk=1,constr_homology
+          if(l_homo(kk,ii)) then 
+            min_odl=distancek(kk)
+            exit
+          endif
+         enddo
+         do kk=1,constr_homology
+          if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) 
+     &              min_odl=distancek(kk)
+         enddo
+
 c        write (iout,* )"min_odl",min_odl
 #ifdef DEBUG
          write (iout,*) "ij dij",i,j,dij

diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F
index a504ff7..de172bd 100644 (file)
--- a/source/unres/src_MD-M/readpdb.F
+++ b/source/unres/src_MD-M/readpdb.F
@@ -1006,6 +1006,7 @@ c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
       endif
       endif
         do j=1,3
+          chomo(j,i,k)=c(j,i)
           cref(j,i,cou)=c(j,i)
           cref(j,i+nres,cou)=c(j,i+nres)
           if (i.le.nres) then

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 7be445d..f944233 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -2580,9 +2580,19 @@ c Alternative: reading from input
       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
- 
+      call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
       call readi(controlcard,"HOMOL_NSET",homol_nset,1)       
       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
+      start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
+      if(.not.read2sigma.and.start_from_model) then
+          write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
+          start_from_model=.false.
+      endif
+      if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
+      if(start_from_model .and. rest) then 
+        write(iout,*) 'START_FROM_MODELS is OFF'
+        write(iout,*) 'remove restart keyword from input'
+      endif
       if (homol_nset.gt.1)then
          call card_concat(controlcard)
          read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
@@ -2679,8 +2689,8 @@ c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
      &                                idomain_tmp
              i_tmp=i_tmp+nnt-1
              idomain(k,i_tmp)=idomain_tmp
-             rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0)
-             rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0)
+             rescore(k,i_tmp)=rescore_tmp
+             rescore2(k,i_tmp)=rescore2_tmp
             else
              idomain(k,irec)=1
              read (ientin,*,end=1401) rescore_tmp
@@ -2697,8 +2707,15 @@ c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
           do i = nnt,nct-2 
             do j=i+2,nct 
 
-            if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
-                  
+              x12=c(1,i)-c(1,j)
+              y12=c(2,i)-c(2,j)
+              z12=c(3,i)-c(3,j)
+              distal=dsqrt(x12*x12+y12*y12+z12*z12) 
+c              write (iout,*) k,i,j,distal,dist2_cut
+
+            if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+     &            .and. distal.le.dist2_cut ) then
+
               ii=ii+1
               ii_in_use(ii)=1
               l_homo(k,ii)=.true.
@@ -2708,13 +2725,6 @@ c             write (iout,*) "i",i," j",j," constr_homology",
 c    &                       constr_homology
               ires_homo(ii)=i
               jres_homo(ii)=j
-c
-c Attempt to replace dist(i,j) by its definition in ...
-c
-              x12=c(1,i)-c(1,j)
-              y12=c(2,i)-c(2,j)
-              z12=c(3,i)-c(3,j)
-              distal=dsqrt(x12*x12+y12*y12+z12*z12)
               odl(k,ii)=distal
               if (read2sigma) then
                 sigma_odl(k,ii)=0
@@ -2726,30 +2736,18 @@ c
      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
               else
                 if (odl(k,ii).le.dist_cut) then
-                 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
-c                sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+                 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) 
                 else
 #ifdef OLDSIGMA
-                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
+                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 #else
-                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
+                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 #endif
                 endif
               endif
-c   Following expr replaced by a positive exp argument
-c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
-c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
-
-c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
-c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
-c
-              sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
-c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
-c
-c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
-c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
+              sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
             else
               ii=ii+1
               l_homo(k,ii)=.false.
@@ -2784,9 +2782,9 @@ c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 
-            sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
-     &                     rescore(k,i-2)+rescore(k,i-3)
-            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
+            sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+     &                     rescore(k,i-2)+rescore(k,i-3))/4.0
+c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
@@ -2821,9 +2819,9 @@ c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
-             sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
-     &                        rescore(k,i-2)
-             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
+             sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
+     &                        rescore(k,i-2))/3.0
+c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 
 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
@@ -2908,24 +2906,24 @@ 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(2f8.2,1x,l1,4x))') 
+       write(iout,'(3i5,100(2f8.2,1x,l1,4x))') 
      &  ii,ires_homo(ii),jres_homo(ii),
      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
      &  ki=1,constr_homology)
        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),
+        write (iout,'(i5,100(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
-        write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+        write (iout,'(i5,100(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
-        write(iout,'(i5,10(4f8.2,4x))') i,
+        write(iout,'(i5,100(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