wham SAXS cutoff

[unres.git] / source / wham / src-M / molread_zs.F

diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F
index 0c1a211..df062e5 100644 (file)
--- a/source/wham/src-M/molread_zs.F
+++ b/source/wham/src-M/molread_zs.F
@@ -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
 
@@ -360,7 +430,7 @@ 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
@@ -463,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)
@@ -688,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