X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Fcluster%2Fwham%2Fsrc-M%2Freadrtns.F;h=60b1ec2d29200f5b5369c9b0098f2d1aafc2a373;hb=c0b48650276909d82678d8b4263049fd1ae77887;hp=2573b962d4773f06675a1a933857d428cb1b646e;hpb=036299742bb15ca3871e0907c6e6e3ea6fe96e0f;p=unres.git

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index 2573b96..60b1ec2 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'
@@ -91,7 +92,8 @@ C long axis of side chain
       min_var=(index(controlcard,'MINVAR').gt.0)
       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
-      call readi(controlcard,'NCUT',ncut,1)
+      call readi(controlcard,'NCUT',ncut,0)
+      call readi(controlcard,'NCLUST',nclust,5)
       call readi(controlcard,'SYM',symetr,1)
       write (iout,*) 'sym', symetr
       call readi(controlcard,'NSTART',nstart,0)
@@ -102,7 +104,8 @@ C long axis of side chain
       lgrp=(index(controlcard,'LGRP').gt.0)
       caonly=(index(controlcard,'CA_ONLY').gt.0)
       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
-      call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
+      if (ncut.gt.0) 
+     & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
       call readi(controlcard,'IOPT',iopt,2) 
       lside = index(controlcard,"SIDE").gt.0
       efree = index(controlcard,"EFREE").gt.0
@@ -125,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
@@ -166,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)
@@ -186,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)
@@ -414,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)
@@ -955,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)
 
@@ -1103,7 +1179,7 @@ c
         tpl_k_rescore="template"//kic2//".sco"
 
         unres_pdb=.false.
-        call readpdb
+        call readpdb_template(k)
         do i=1,2*nres
           do j=1,3
             crefjlee(j,i)=c(j,i)