SAXS cutoff correct gradiend and cutoff scaling

[unres.git] / source / unres / src_MD-M / readrtns_CSA.F

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 2cbb825..70bb647 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -98,6 +98,13 @@ c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 C Set up the time limit (caution! The time must be input in minutes!)
       read_cart=index(controlcard,'READ_CART').gt.0
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+      write (iout,*) "constr_dist",constr_dist
+      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
       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
       call readi(controlcard,'SYM',symetr,1)
       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
@@ -619,6 +626,7 @@ C Read weights of the subsequent energy terms.
        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
        call reada(weightcard,'TEMP0',temp0,300.0d0)
        call reada(weightcard,'WLT',wliptran,0.0D0)
+       call reada(weightcard,'WSAXS',wsaxs,1.0D0)
        if (index(weightcard,'SOFT').gt.0) ipot=6
 C 12/1/95 Added weight for the multi-body term WCORR
        call reada(weightcard,'WCORRH',wcorr,1.0D0)
@@ -642,6 +650,7 @@ C 12/1/95 Added weight for the multi-body term WCORR
        weights(17)=wbond
        weights(18)=scal14
        weights(21)=wsccor
+       weights(25)=wsaxs
       if(me.eq.king.or..not.out1file)
      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
@@ -803,12 +812,16 @@ c        print *,'Finished reading pdb data'
               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
               nsi=nsi+1
             enddo
+c AL 7/10/16
+c Calculalte only the coordinates of the current sidechain; no need to rebuild
+c whole chain
+            call locate_side_chain(i)
             if(fail) write(iout,*)'Adding sidechain failed for res ',
      &              i,' after ',nsi,' trials'
           endif
          enddo
 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
-         call chainbuild
+c         call chainbuild
         endif  
       endif
       if (indpdb.eq.0) then
@@ -869,6 +882,7 @@ C 8/13/98 Set limits to generating the dihedral angles
         phibound(2,i)=pi
       enddo
       read (inp,*) ndih_constr
+      write (iout,*) "ndish_constr",ndih_constr
       if (ndih_constr.gt.0) then
         read (inp,*) ftors
         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
@@ -998,6 +1012,8 @@ c        write (iout,*) "After read_dist_constr nhpb",nhpb
         enddo
         endif
       endif
+      write (iout,*) "calling read_saxs_consrtr",nsaxs
+      if (nsaxs.gt.0) call read_saxs_constr
 
 
       if (constr_homology.gt.0) then
@@ -2120,6 +2136,8 @@ C Get parameter filenames and open the parameter files.
       open (ielep,file=elename,status='old')
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old')
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
      &  readonly)
@@ -2447,6 +2465,67 @@ CCCC NOW PROPERTIES FOR AFM
       end
 
 c-------------------------------------------------------------------------------
+      subroutine read_saxs_constr
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      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_dist_constr
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'