cluster wham SACS cutoff

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

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index a3feaa5..60b1ec2 100644 (file)
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -13,10 +13,11 @@ C
       include 'COMMON.CHAIN'
       include 'COMMON.HEADER'
       include 'COMMON.FFIELD'
-      include 'COMMON.FREE'
       include 'COMMON.INTERACT'
       include "COMMON.SPLITELE"
       include 'COMMON.SHIELD'
+      include 'COMMON.FREE'
+      include 'COMMON.LANGEVIN'
       character*320 controlcard,ucase
 #ifdef MPL
       include 'COMMON.INFO'
@@ -127,6 +128,12 @@ C long axis of side chain
       print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
       print_homology_models=
      & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
+      call readi(controlcard,'NSAXS',nsaxs,0)
+      call readi(controlcard,'SAXS_MODE',saxs_mode,0)
+      call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
+      call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
+      write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+     &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
       if (min_var) iopt=1
       return
       end
@@ -168,8 +175,7 @@ C
 C Read weights of the subsequent energy terms.
       call card_concat(weightcard)
       write(iout,*) weightcard
-C      call reada(weightcard,'WSC',wsc,1.0d0)
-      write(iout,*) wsc
+      call reada(weightcard,'WSC',wsc,1.0d0)
       call reada(weightcard,'WLONG',wsc,wsc)
       call reada(weightcard,'WSCP',wscp,1.0d0)
       call reada(weightcard,'WELEC',welec,1.0D0)
@@ -188,6 +194,7 @@ C      call reada(weightcard,'WSC',wsc,1.0d0)
       call reada(weightcard,'WTORD',wtor_d,1.0D0)
       call reada(weightcard,'WANG',wang,1.0D0)
       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+      call reada(weightcard,'WSAXS',wsaxs,0.0D0)
       call reada(weightcard,'SCAL14',scal14,0.4D0)
       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
@@ -416,6 +423,8 @@ C        enddo
       if (nstart.lt.nnt) nstart=nnt
       if (nend.gt.nct .or. nend.eq.0) nend=nct
       write (iout,*) "nstart",nstart," nend",nend
+      write (iout,*) "calling read_saxs_consrtr",nsaxs
+      if (nsaxs.gt.0) call read_saxs_constr
       nres0=nres
       if (constr_homology.gt.0) then
         call read_constr_homology(print_homology_restraints)
@@ -957,6 +966,71 @@ C      endif
       return
       end
 
+c-------------------------------------------------------------------------------
+      subroutine read_saxs_constr
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#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(lprn)