X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Freadrtns.F;h=e5c31cf33efc9a894ed18efcc2407c6e125bfe2a;hb=17f6ba1f1fcc75e0823f1cfb16fa99c11f606d6d;hp=a3feaa570b5679592a1caa0be11f8db150f27de4;hpb=ba8f40a2336651a42eb42e9d3379d62aa11bb5d7;p=unres.git

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index a3feaa5..e5c31cf 100644
--- 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,11 @@ 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)
+      write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+     &   SAXS_MODE," SCAL_RAD",scal_rad
       if (min_var) iopt=1
       return
       end
@@ -168,8 +174,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 +193,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 +422,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 +965,66 @@ 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
+      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)