unres_package_Oct_2016 from emilial
[unres4.git] / source / wham / wham.f90
diff --git a/source/wham/wham.f90 b/source/wham/wham.f90
new file mode 100644 (file)
index 0000000..fcf1d15
--- /dev/null
@@ -0,0 +1,372 @@
+      program wham_multparm
+!      program WHAM_multparm
+! Creation/update of the database of conformations
+      use wham_data
+      use io_wham
+      use io_database
+      use wham_calc
+      use ene_calc
+      use conform_compar
+      use work_part
+!
+      use io_units
+      use control_data, only:indpdb
+#ifdef MPI
+      use mpi_data
+!      use mpi_
+#endif
+      use control, only:initialize
+!el      use io_config, only:parmread
+!
+#ifndef ISNAN
+      external proc_proc
+#ifdef WINPGI
+!MS$ATTRIBUTES C ::  proc_proc
+#endif
+#endif
+!
+!el#ifndef ISNAN
+!el      external proc_proc
+!el#endif
+!el#ifdef WINPGI
+!elcMS$ATTRIBUTES C ::  proc_proc
+!el#endif
+!      include "DIMENSIONS"
+!      include "DIMENSIONS.ZSCOPT"
+!      include "DIMENSIONS.FREE"
+!      implicit none
+#ifdef MPI
+!      include "COMMON.MPI"
+!      use mpi_data
+      include "mpif.h"
+      integer :: IERROR,ERRCODE
+#endif
+!      include "COMMON.IOUNITS"
+!      include "COMMON.FREE"
+!      include "COMMON.CONTROL"
+!      include "COMMON.ALLPARM"
+!      include "COMMON.PROT"
+      real(kind=8) :: rr !,x(max_paropt)
+      integer :: idumm
+      integer :: i,ipar,islice
+
+!el      run_wham=.true.
+!#define WHAM_RUN
+!      call alloc_wham_arrays
+!write(iout,*) "after alloc wham"
+#ifdef MPI
+      call MPI_Init( IERROR )
+      call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
+      call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
+      Master = 0
+      if (ierror.gt.0) then
+        write(iout,*) "SEVERE ERROR - Can't initialize MPI."
+        call mpi_finalize(ierror)
+        stop
+      endif
+!el      if (nprocs.gt.MaxProcs+1) then
+!el        write (2,*) "Error - too many processors",&
+!el         nprocs,MaxProcs+1
+!el        write (2,*) "Increase MaxProcs and recompile"
+!el        call MPI_Finalize(IERROR)
+!el        stop
+!el      endif
+#endif
+! NaNQ initialization
+#ifndef ISNAN
+      i=-1
+      rr=dacos(100.0d0)
+#ifdef WINPGI
+      idumm=proc_proc(rr,i)
+#else
+      call proc_proc(rr,i)
+#endif
+#endif
+!write(iout,*) "before init"
+      call initialize
+!write(iout,*)"after init"
+      call openunits
+!write(iout,*)"after open ui"
+      call cinfo
+!write(iout,*)"after cinfo"
+      call read_general_data(*10)
+!write(iout,*)"after read_gen"
+      call flush(iout)
+      call molread(*10)
+!write(iout,*)"after molread"
+      call flush(iout)
+#ifdef MPI 
+      write (iout,*) "Calling proc_groups"
+      call proc_groups
+      write (iout,*) "proc_groups exited"
+      call flush(iout)
+#endif
+!el----------
+      call alloc_wham_arrays
+!el----------
+      do ipar=1,nParmSet
+        write (iout,*) "Calling parmread",ipar
+        call parmread(ipar,*10)
+        if (.not.separate_parset) then
+          call store_parm(ipar)
+          write (iout,*) "Finished storing parameters",ipar
+        else if (ipar.eq.myparm) then
+          call store_parm(1)
+          write (iout,*) "Finished storing parameters",ipar
+        endif
+        call flush(iout)
+      enddo
+      call read_efree(*10)
+      write (iout,*) "Finished READ_EFREE"
+      call flush(iout)
+      call read_protein_data(*10)
+      write (iout,*) "Finished READ_PROTEIN_DATA"
+      call flush(iout)
+      if (indpdb.gt.0) then
+        call promienie
+        call read_compar
+        call read_ref_structure(*10)
+!write(iout,*)"before proc_cont, define frag"
+        call proc_cont
+        call fragment_list
+        if (constr_dist.gt.0) call read_dist_constr
+      endif
+      write (iout,*) "Begin read_database"
+      call flush(iout)
+      call read_database(*10)
+      write (iout,*) "Finished read_database"
+      call flush(iout)
+      if (separate_parset) nparmset=1
+      do islice=1,nslice
+        if (ntot(islice).gt.0) then
+#ifdef MPI 
+          call work_partition(islice,.true.)
+          write (iout,*) "work_partition OK"
+          call flush(iout)
+#endif
+          write (iout,*) "call enecalc",islice,nslice
+          call enecalc(islice,*10)
+          write (iout,*) "enecalc OK"
+          call flush(iout)
+          call WHAMCALC(islice,*10)
+          write (iout,*) "wham_calc OK"
+          call flush(iout)
+          call write_dbase(islice,*10)
+          write (iout,*) "write_dbase OK"
+          call flush(iout)
+          if (ensembles.gt.0) then
+            call make_ensembles(islice,*10)
+            write (iout,*) "make_ensembles OK"
+            call flush(iout)
+          endif
+        endif
+      enddo
+#ifdef MPI
+      call MPI_Finalize( IERROR )
+#endif
+      stop
+   10 write (iout,*) "Error termination of the program"
+#ifdef MPI
+      call MPI_Finalize( IERROR )
+#endif
+      stop
+      end program wham_multparm
+!------------------------------------------------------------------------------
+!
+!------------------------------------------------------------------------------
+#ifdef MPI
+      subroutine proc_groups
+! Split the processors into the Master and Workers group, if needed.
+      use io_units
+      use MPI_data
+      use wham_data
+
+      implicit none
+!      include "DIMENSIONS"
+!      include "DIMENSIONS.ZSCOPT"
+!      include "DIMENSIONS.FREE"
+!      include "COMMON.IOUNITS"
+!      include "COMMON.MPI"
+!      include "COMMON.FREE"
+      include "mpif.h"
+      integer :: n,chunk,i,j,ii,remainder
+      integer :: kolorW,key,ierror,errcode
+      logical :: lprint
+      lprint=.true.
+!
+! Split the communicator if independent runs for different parameter
+! sets will be performed.
+!
+      if (nparmset.eq.1 .or. .not.separate_parset) then
+        WHAM_COMM = MPI_COMM_WORLD
+      else if (separate_parset) then
+        if (nprocs.lt.nparmset) then
+          write (iout,*) &
+       "*** Cannot split parameter sets for fewer processors than sets",&
+        nprocs,nparmset
+          call MPI_Finalize(ierror)
+          stop
+        endif 
+        write (iout,*) "nparmset",nparmset
+        nprocs = nprocs/nparmset
+        kolorW = me/nprocs
+        key = mod(me,nprocs)
+        write (iout,*) "My old rank",me," kolor",kolorW," key",key
+        call MPI_Comm_split(MPI_COMM_WORLD,kolorW,key,WHAM_COMM,ierror)
+        call MPI_Comm_size(WHAM_COMM,nprocs,ierror)
+        call MPI_Comm_rank(WHAM_COMM,me,ierror)
+        write (iout,*) "My new rank",me," comm size",nprocs
+        write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
+         " WHAM_COMM",WHAM_COMM
+        myparm=kolorW+1
+        write (iout,*) "My parameter set is",myparm
+        call flush(iout)
+      else
+        myparm=nparmset
+      endif
+      Me1 = Me
+      Nprocs1 = Nprocs
+      return
+      end subroutine proc_groups
+#endif
+!------------------------------------------------------------------------------
+#ifdef AIX
+      subroutine flush(iu)
+      call flush_(iu)
+      return
+      end subroutine flush
+#endif
+!-----------------------------------------------------------------------------
+      subroutine promienie(*)
+
+      use io_units
+      use names
+      use io_base, only:ucase
+      use energy_data, only:sigma0,dsc,dsc_inv
+      use wham_data
+      use w_compar_data
+      implicit none
+!      include 'DIMENSIONS'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CONTPAR'
+!      include 'COMMON.LOCAL'
+      integer ::i,j
+      real(kind=8) :: facont=1.569D0  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
+      character(len=8) :: contfunc
+      character(len=8) :: contfuncid(5)=reshape((/'GB      ',&
+                      'DIST    ','CEN     ','ODC     ','SIG     '/),shape(contfuncid))
+!el      character(len=8) ucase
+      call getenv("CONTFUNC",contfunc)
+      contfunc=ucase(contfunc)
+      do icomparfunc=1,5
+        if (contfunc.eq.contfuncid(icomparfunc)) goto 10
+      enddo     
+   10 continue
+      write (iout,*) "Sidechain contact function is ",contfunc,&
+        "icomparfunc",icomparfunc 
+      do i=1,ntyp
+        do j=1,ntyp
+          if (icomparfunc.lt.3) then
+            read(isidep1,*) chi_comp(i,j),chip_comp(i,j),sig_comp(i,j),&
+             sc_cutoff(i,j)
+          else if (icomparfunc.lt.5) then
+            read(isidep1,*) sc_cutoff(i,j)
+          else if (icomparfunc.eq.5) then
+            sc_cutoff(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)*facont
+          else
+            write (iout,*) "Error - Unknown contact function"
+            return 1
+          endif
+        enddo
+      enddo
+      close (isidep1)
+      do i=1,ntyp1
+        if (i.eq.10 .or. i.eq.ntyp1) then
+          dsc_inv(i)=0.0d0
+        else
+          dsc_inv(i)=1.0d0/dsc(i)
+        endif
+      enddo
+      return
+      end subroutine promienie
+!-----------------------------------------------------------------------------
+      subroutine alloc_wham_arrays
+
+      use names
+      use geometry_data, only:nres
+      use energy_data, only:maxcont
+      use wham_data
+      use w_compar_data
+      integer :: i,j,k,l
+!-------------------------
+! COMMON.FREE
+!      common /wham/
+      allocate(stot(nslice)) !(maxslice)
+      do i=1,nslice
+        stot(i)=0
+      enddo
+      allocate(Kh(nQ,MaxR,MaxT_h,nParmSet),q0(nQ,MaxR,MaxT_h,nParmSet))!(MaxQ,MaxR,MaxT_h,max_parm)
+      allocate(f(maxR,maxT_h,nParmSet)) !(maxR,maxT_h,max_parm)
+      allocate(beta_h(maxT_h,nParmSet)) !(MaxT_h,max_parm)
+      allocate(nR(maxT_h,nParmSet),nRR(maxT_h,nParmSet)) !(maxT_h,max_parm)
+      allocate(snk(MaxR,MaxT_h,nParmSet,nSlice)) !(MaxR,MaxT_h,max_parm,MaxSlice)
+!      do i=1,MaxR
+!        do j=1,MaxT_h
+!          do k=1,nParmSet
+!            do l=1,nSlice
+!              snk(i,j,k,l)=0
+!            enddo
+!          enddo
+!        enddo
+!      enddo
+
+      allocate(totraj(maxR,nParmSet)) !(maxR,max_parm)
+
+      allocate(nT_h(nParmSet))!(max_parm)
+      allocate(replica(nParmSet))
+      allocate(umbrella(nParmSet))
+      allocate(read_iset(nParmSet))
+!      allocate(nT_h(nParmSet))
+!-------------------------
+! COMMON.PROT
+!      common /protein/
+      allocate(ntot(nslice))  !(maxslice)
+!      allocatable :: isampl        !(max_parm)
+!-------------------------
+! COMMON.PROTFILES
+!      common /protfil/
+      allocate(protfiles(maxfile_prot,2,MaxR,MaxT_h,nParmSet)) !(maxfile_prot,2,MaxR,MaxT_h,Max_Parm)
+      allocate(nfile_bin(MaxR,MaxT_h,nParmSet))
+      allocate(nfile_asc(MaxR,MaxT_h,nParmSet))
+      allocate(nfile_cx(MaxR,MaxT_h,nParmSet))
+      allocate(rec_start(MaxR,MaxT_h,nParmSet))
+      allocate(rec_end(MaxR,MaxT_h,nParmSet)) !(MaxR,MaxT_h,Max_Parm)
+!-------------------------
+! COMMON.OBCINKA
+!      common /obcinka/
+      allocate(time_start_collect(maxR,MaxT_h,nParmSet))
+      allocate(time_end_collect(maxR,MaxT_h,nParmSet)) !(maxR,MaxT_h,Max_Parm)
+!-------------------------
+! COMMON.CONTPAR
+!      common /contpar/
+      allocate(sig_comp(ntyp,ntyp),chi_comp(ntyp,ntyp),&
+        chip_comp(ntyp,ntyp),sc_cutoff(ntyp,ntyp)) !(ntyp,ntyp)
+!-------------------------
+! COMMON.PEPTCONT
+!      common /peptcont/
+      allocate(icont_pept_ref(2,maxcont)) !(2,maxcont)
+!      allocate(ncont_frag_ref()) !(mmaxfrag)
+!      allocate(icont_frag_ref(2,maxcont)) !(2,maxcont,mmaxfrag)
+      allocate(isec_ref(nres)) !(maxres)
+!-------------------------
+! COMMON.VAR
+! Angles from experimental structure
+!      common /varref/
+      allocate(vbld_ref(nres),theta_ref(nres),&
+        phi_ref(nres),alph_ref(nres),omeg_ref(nres)) !(maxres)
+!-------------------------
+      end subroutine alloc_wham_arrays
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------