wham SAXS cutoff
[unres.git] / source / wham / src-M / molread_zs.F
index dee1afc..df062e5 100644 (file)
@@ -137,6 +137,9 @@ C        enddo
       if (itype(1).eq.ntyp1) nnt=2
       if (itype(nres).eq.ntyp1) nct=nct-1
       write(iout,*) 'NNT=',NNT,' NCT=',NCT
+      write (iout,*) "calling read_saxs_consrtr",nsaxs
+      if (nsaxs.gt.0) call read_saxs_constr
+
       if (constr_homology.gt.0) then
 c       write (iout,*) "About to call read_constr_homology"
 c       call flush(iout)
@@ -324,6 +327,73 @@ c      print *,"energy",energ," iscor",iscor
       return
    10 return1
       end
+c-------------------------------------------------------------------------------
+      subroutine read_saxs_constr
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SBRIDGE'
+      double precision cm(3)
+c      read(inp,*) nsaxs
+      write (iout,*) "Calling read_saxs nsaxs",nsaxs
+      call flush(iout)
+      if (saxs_mode.eq.0) then
+c SAXS distance distribution
+      do i=1,nsaxs
+        read(inp,*) distsaxs(i),Psaxs(i)
+      enddo
+      Cnorm = 0.0d0
+      do i=1,nsaxs
+        Cnorm = Cnorm + Psaxs(i)
+      enddo
+      write (iout,*) "Cnorm",Cnorm
+      do i=1,nsaxs
+        Psaxs(i)=Psaxs(i)/Cnorm
+      enddo
+      write (iout,*) "Normalized distance distribution from SAXS"
+      do i=1,nsaxs
+        write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
+      enddo
+      Wsaxs0=0.0d0
+      do i=1,nsaxs
+        Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
+      enddo
+      write (iout,*) "Wsaxs0",Wsaxs0
+      else
+c SAXS "spheres".
+      do i=1,nsaxs
+        read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
+      enddo
+      do j=1,3
+        cm(j)=0.0d0
+      enddo
+      do i=1,nsaxs
+        do j=1,3
+          cm(j)=cm(j)+Csaxs(j,i)
+        enddo
+      enddo
+      do j=1,3
+        cm(j)=cm(j)/nsaxs
+      enddo
+      do i=1,nsaxs
+        do j=1,3
+          Csaxs(j,i)=Csaxs(j,i)-cm(j)
+        enddo
+      enddo
+      write (iout,*) "SAXS sphere coordinates"
+      do i=1,nsaxs
+        write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
+      enddo
+      endif
+      return
+      end
 c====-------------------------------------------------------------------
       subroutine read_constr_homology
 
@@ -339,6 +409,7 @@ c====-------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
+      include 'include_unres/COMMON.NAMES'
       include 'COMMON.HOMRESTR'
 c
 c For new homol impl
@@ -359,13 +430,15 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       logical lprn /.true./
       integer ilen
       external ilen
-      logical unres_pdb
+      logical unres_pdb,liiflag
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
-      double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+      double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
+     &                 rescore3_tmp
       double precision, dimension (max_template,maxres) :: rescore
       double precision, dimension (max_template,maxres) :: rescore2
+      double precision, dimension (max_template,maxres) :: rescore3
       character*24 tpl_k_rescore
 c -----------------------------------------------------------------
 c Reading multiple PDB ref structures and calculation of retraints
@@ -410,8 +483,6 @@ cd      call flush(iout)
 c
 c  New
 c
-      lim_theta=0
-      lim_xx=0
 c
 c  Reading HM global scores (prob not required)
 c
@@ -462,7 +533,7 @@ c
         tpl_k_rescore="template"//kic2//".sco"
 
         unres_pdb=.false.
-        call readpdb
+        call readpdb_template(k)
 cref        do i=1,2*nres
 cref          do j=1,3
 cref            crefjlee(j,i)=c(j,i)
@@ -494,14 +565,18 @@ c
 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
           open (ientin,file=tpl_k_rescore,status='old')
           if (nnt.gt.1) rescore(k,1)=0.0d0
-          do irec=nnt,maxdim ! loop for reading res sim 
+          do irec=nnt,nct ! loop for reading res sim 
             if (read2sigma) then
              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
-     &                                idomain_tmp
+     &                                rescore3_tmp,idomain_tmp
              i_tmp=i_tmp+nnt-1
              idomain(k,i_tmp)=idomain_tmp
              rescore(k,i_tmp)=rescore_tmp
              rescore2(k,i_tmp)=rescore2_tmp
+             rescore3(k,i_tmp)=rescore3_tmp
+             write(iout,'(a7,i5,3f10.5,i5)') "rescore",
+     &                      i_tmp,rescore2_tmp,rescore_tmp,
+     &                                rescore3_tmp,idomain_tmp
             else
              idomain(k,irec)=1
              read (ientin,*,end=1401) rescore_tmp
@@ -601,7 +676,8 @@ 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)*
 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 c   Instead of res sim other local measure of b/b str reliability possible
-            sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+            if (sigma_dih(k,i).ne.0)
+     &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
           enddo
           lim_dih=nct-nnt-2 
@@ -633,14 +709,14 @@ 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
 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))
+             if (sigma_theta(k,i).ne.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)*
 c                             rescore(k,i-2) !  right expression ?
 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
           enddo
         endif
-        lim_theta=nct-nnt-1 
 
         if (waga_d.gt.0.0d0) then
 c       open (ientin,file=tpl_k_sigma_d,status='old')
@@ -666,15 +742,14 @@ c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
-               sigma_d(k,i)=rescore(k,i) !  right expression ?
-               sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+               sigma_d(k,i)=rescore3(k,i) !  right expression ?
+               if (sigma_d(k,i).ne.0)
+     &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 
 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?
           enddo
-          lim_xx=nct-nnt+1 
         endif
       enddo
 c
@@ -683,22 +758,34 @@ c shift data in all arrays
 c
       if (waga_dist.ne.0.0d0) then
         ii=0
+        liiflag=.true.
         do i=nnt,nct-2 
          do j=i+2,nct 
           ii=ii+1
-          if (ii_in_use(ii).eq.0) then 
-             do ki=ii,lim_odl-1
-              ires_homo(ki)=ires_homo(ki+1)
-              jres_homo(ki)=jres_homo(ki+1)
-              ii_in_use(ki)=ii_in_use(ki+1)
+          if (ii_in_use(ii).eq.0.and.liiflag) then
+            liiflag=.false.
+            iistart=ii
+          endif
+          if (ii_in_use(ii).ne.0.and..not.liiflag.or.
+     &                   .not.liiflag.and.ii.eq.lim_odl) then
+             if (ii.eq.lim_odl) then
+              iishift=ii-iistart+1
+             else
+              iishift=ii-iistart
+             endif
+             liiflag=.true.
+             do ki=iistart,lim_odl-iishift
+              ires_homo(ki)=ires_homo(ki+iishift)
+              jres_homo(ki)=jres_homo(ki+iishift)
+              ii_in_use(ki)=ii_in_use(ki+iishift)
               do k=1,constr_homology
-               odl(k,ki)=odl(k,ki+1)
-               sigma_odl(k,ki)=sigma_odl(k,ki+1)
-               l_homo(k,ki)=l_homo(k,ki+1)
+               odl(k,ki)=odl(k,ki+iishift)
+               sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
+               l_homo(k,ki)=l_homo(k,ki+iishift)
               enddo
              enddo
-             ii=ii-1
-             lim_odl=lim_odl-1
+             ii=ii-iishift
+             lim_odl=lim_odl-iishift
           endif
          enddo
         enddo
@@ -723,17 +810,19 @@ cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
      &  ki=1,constr_homology)
        enddo
        write (iout,*) "Dihedral angle restraints from templates"
-       do i=nnt+3,lim_dih
-        write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+       do i=nnt+3,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+       do i=nnt+2,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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
+       do i=nnt,nct
         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)