restoring read2sigma code after wrong merge
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 31 Jul 2015 19:33:28 +0000 (21:33 +0200)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 31 Jul 2015 19:33:28 +0000 (21:33 +0200)
source/unres/src_MD/readrtns.F

index dd5536b..ecbb19f 100644 (file)
@@ -2726,8 +2726,9 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
-      double precision rescore_tmp,x12,y12,z12
+      double precision rescore_tmp,x12,y12,z12,rescore2_tmp
       double precision, dimension (max_template,maxres) :: rescore
+      double precision, dimension (max_template,maxres) :: rescore2
       character*24 pdbfile,tpl_k_rescore
 c -----------------------------------------------------------------
 c Reading multiple PDB ref structures and calculation of retraints
@@ -2745,6 +2746,7 @@ c Alternative: reading from input
       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
  
       call readi(controlcard,"HOMOL_NSET",homol_nset,1)       
+      read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
       if (homol_nset.gt.1)then
          call card_concat(controlcard)
          read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
         write(kic2,'(bz,i2.2)') k
 
         tpl_k_rescore="template"//kic2//".sco"
-c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
-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
+        if (read2sigma) then
+          call readpdb_template(k)
+        else
+          call readpdb
+        endif
 c
 c     Distance restraints
 c
@@ -2853,24 +2855,11 @@ c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
             endif
           enddo  
  1401   continue
-          close (ientin)        
-c         open (ientin,file=tpl_k_sigma_odl,status='old')
-c         do irec=1,maxdim ! loop for reading sigma_odl
-c            read (ientin,*,end=1401) i, j, 
-c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
-c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
-c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k) 
-c         enddo
-c 1401   continue
-c         close (ientin)
+        close (ientin)        
         if (waga_dist.ne.0.0d0) then
           ii=0
-          do i = nnt,nct-2 ! right? without parallel.
-            do j=i+2,nct ! right?
-c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology 
-c           do j=i+2,nres ! ibid
-c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
-c           do j=i+2,nct ! ibid
+          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
                   
               z12=c(3,i)-c(3,j)
               distal=dsqrt(x12*x12+y12*y12+z12*z12)
               odl(k,ii)=distal
-c
-c             odl(k,ii)=dist(i,j)
-c             write (iout,*) "dist(",i,j,") =",dist(i,j)
-c             write (iout,*) "distal = ",distal
-c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
-c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
-c    &                      "rescore(",k,j,") =",rescore(k,j)
-c
-c  Calculation of sigma from res sim
-c
-c             if (odl(k,ii).le.6.0d0) then
-c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
-c    Other functional forms possible depending on odl(k,ii), eg.
-c
-            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)
-            else
+              if (read2sigma) then
+                sigma_odl(k,ii)=0
+                do ik=i,j
+                 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+                enddo
+                sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+                if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
+     &        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)
+                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))* ! sigma ~ rescore ~ error 
      &                      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))* ! sigma ~ rescore ~ error 
      &                      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)
-            endif
 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)
@@ -2935,10 +2919,7 @@ c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of us
               l_homo(k,ii)=.false.
             endif
             enddo
-c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
           enddo
-c         lim_odl=ii
-c         if (constr_homology.gt.0) call homology_partition
         lim_odl=ii
         endif
 c
@@ -2954,10 +2935,11 @@ c    &                            sigma_dih(k,i+nnt-1)
 c         enddo
 c1402   continue
 c         close (ientin)
-          do i = nnt+3,nct ! right? without parallel.
-c         do i=1,nres ! alternative for bounds acc to readpdb?
-c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
+          do i = nnt+3,nct 
+            if (idomain(k,i).eq.0) then 
+               sigma_dih(k,i)=0.0
+               cycle
+            endif
             dih(k,i)=phiref(i) ! right?
 c           read (ientin,*) sigma_dih(k,i) ! original variant
 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
@@ -2993,6 +2975,10 @@ c         close (ientin)
           do i = nnt+2,nct ! right? without parallel.
 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 c         do i=ithet_start,ithet_end ! with FG parallel.
+             if (idomain(k,i).eq.0) then  
+              sigma_theta(k,i)=0.0
+              cycle
+             endif
              thetatpl(k,i)=thetaref(i)
 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
@@ -3000,7 +2986,7 @@ 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)/3.0
+     &                        rescore(k,i-2))/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)*
@@ -3018,12 +3004,15 @@ c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use
 c    &                          sigma_d(k,i+nnt-1)
 c         enddo
 c1404   continue
-          close (ientin)
 
           do i = nnt,nct ! right? without parallel.
 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 c         do i=loc_start,loc_end ! with FG parallel.
-             if (itype(i).eq.10) goto 1 ! right?
+               if (itype(i).eq.10) cycle 
+               if (idomain(k,i).eq.0 ) then 
+                  sigma_d(k,i)=0.0
+                  cycle
+               endif
                xxtpl(k,i)=xxref(i)
                yytpl(k,i)=yyref(i)
                zztpl(k,i)=zzref(i)
@@ -3038,11 +3027,9 @@ c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 c              read (ientin,*) sigma_d(k,i) ! 1st variant
                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
-    1     continue
           enddo
           lim_xx=nct-nnt+1 
         endif
-        close(ientin)
       enddo
 c
 c remove distance restraints not used in any model from the list