source/wham/wham.f90

   1       program wham_multparm
   2 !      program WHAM_multparm
   3 ! Creation/update of the database of conformations
   4       use wham_data
   5       use io_wham
   6       use io_database
   7       use wham_calc
   8       use ene_calc
   9       use conform_compar
  10       use work_part
  11 !
  12       use io_units
  13       use control_data, only:indpdb
  14 #ifdef MPI
  15       use mpi_data
  16 !      use mpi_
  17 #endif
  18       use control, only:initialize
  19 !el      use io_config, only:parmread
  20 !
  21 #ifndef ISNAN
  22       external proc_proc
  23 #ifdef WINPGI
  24 !MS$ATTRIBUTES C ::  proc_proc
  25 #endif
  26 #endif
  27 !
  28 !el#ifndef ISNAN
  29 !el      external proc_proc
  30 !el#endif
  31 !el#ifdef WINPGI
  32 !elcMS$ATTRIBUTES C ::  proc_proc
  33 !el#endif
  34 !      include "DIMENSIONS"
  35 !      include "DIMENSIONS.ZSCOPT"
  36 !      include "DIMENSIONS.FREE"
  37 !      implicit none
  38 #ifdef MPI
  39 !      include "COMMON.MPI"
  40 !      use mpi_data
  41       include "mpif.h"
  42       integer :: IERROR,ERRCODE
  43 #endif
  44 !      include "COMMON.IOUNITS"
  45 !      include "COMMON.FREE"
  46 !      include "COMMON.CONTROL"
  47 !      include "COMMON.ALLPARM"
  48 !      include "COMMON.PROT"
  49       real(kind=8) :: rr !,x(max_paropt)
  50       integer :: idumm
  51       integer :: i,ipar,islice
  52
  53 !el      run_wham=.true.
  54 !#define WHAM_RUN
  55 !      call alloc_wham_arrays
  56 !write(iout,*) "after alloc wham"
  57 #ifdef MPI
  58       call MPI_Init( IERROR )
  59       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
  60       call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
  61       Master = 0
  62       if (ierror.gt.0) then
  63         write(iout,*) "SEVERE ERROR - Can't initialize MPI."
  64         call mpi_finalize(ierror)
  65         stop
  66       endif
  67 !el      if (nprocs.gt.MaxProcs+1) then
  68 !el        write (2,*) "Error - too many processors",&
  69 !el         nprocs,MaxProcs+1
  70 !el        write (2,*) "Increase MaxProcs and recompile"
  71 !el        call MPI_Finalize(IERROR)
  72 !el        stop
  73 !el      endif
  74 #endif
  75 ! NaNQ initialization
  76 #ifndef ISNAN
  77       i=-1
  78       rr=dacos(100.0d0)
  79 #ifdef WINPGI
  80       idumm=proc_proc(rr,i)
  81 #else
  82       call proc_proc(rr,i)
  83 #endif
  84 #endif
  85 !write(iout,*) "before init"
  86       call initialize
  87 !write(iout,*)"after init"
  88       call openunits
  89 !write(iout,*)"after open ui"
  90       call cinfo
  91 !write(iout,*)"after cinfo"
  92       call read_general_data(*10)
  93 !write(iout,*)"after read_gen"
  94       call flush(iout)
  95       call molread(*10)
  96 !write(iout,*)"after molread"
  97       call flush(iout)
  98 #ifdef MPI
  99       write (iout,*) "Calling proc_groups"
 100       call proc_groups
 101       write (iout,*) "proc_groups exited"
 102       call flush(iout)
 103 #endif
 104 !el----------
 105       call alloc_wham_arrays
 106 !el----------
 107       do ipar=1,nParmSet
 108         write (iout,*) "Calling parmread",ipar
 109         call parmread(ipar,*10)
 110         if (.not.separate_parset) then
 111           call store_parm(ipar)
 112           write (iout,*) "Finished storing parameters",ipar
 113         else if (ipar.eq.myparm) then
 114           call store_parm(1)
 115           write (iout,*) "Finished storing parameters",ipar
 116         endif
 117         call flush(iout)
 118       enddo
 119       call read_efree(*10)
 120       write (iout,*) "Finished READ_EFREE"
 121       call flush(iout)
 122       call read_protein_data(*10)
 123       write (iout,*) "Finished READ_PROTEIN_DATA"
 124       call flush(iout)
 125       if (indpdb.gt.0) then
 126         call promienie
 127         call read_compar
 128         call read_ref_structure(*10)
 129 !write(iout,*)"before proc_cont, define frag"
 130         call proc_cont
 131         call fragment_list
 132         if (constr_dist.gt.0) call read_dist_constr
 133       endif
 134       write (iout,*) "Begin read_database"
 135       call flush(iout)
 136       call read_database(*10)
 137       write (iout,*) "Finished read_database"
 138       call flush(iout)
 139       if (separate_parset) nparmset=1
 140       do islice=1,nslice
 141         if (ntot(islice).gt.0) then
 142 #ifdef MPI
 143           call work_partition(islice,.true.)
 144           write (iout,*) "work_partition OK"
 145           call flush(iout)
 146 #endif
 147           write (iout,*) "call enecalc",islice,nslice
 148           call enecalc(islice,*10)
 149           write (iout,*) "enecalc OK"
 150           call flush(iout)
 151           call WHAMCALC(islice,*10)
 152           write (iout,*) "wham_calc OK"
 153           call flush(iout)
 154           call write_dbase(islice,*10)
 155           write (iout,*) "write_dbase OK"
 156           call flush(iout)
 157           if (ensembles.gt.0) then
 158             call make_ensembles(islice,*10)
 159             write (iout,*) "make_ensembles OK"
 160             call flush(iout)
 161           endif
 162         endif
 163       enddo
 164 #ifdef MPI
 165       call MPI_Finalize( IERROR )
 166 #endif
 167       stop
 168    10 write (iout,*) "Error termination of the program"
 169 #ifdef MPI
 170       call MPI_Finalize( IERROR )
 171 #endif
 172       stop
 173       end program wham_multparm
 174 !------------------------------------------------------------------------------
 175 !
 176 !------------------------------------------------------------------------------
 177 #ifdef MPI
 178       subroutine proc_groups
 179 ! Split the processors into the Master and Workers group, if needed.
 180       use io_units
 181       use MPI_data
 182       use wham_data
 183
 184       implicit none
 185 !      include "DIMENSIONS"
 186 !      include "DIMENSIONS.ZSCOPT"
 187 !      include "DIMENSIONS.FREE"
 188 !      include "COMMON.IOUNITS"
 189 !      include "COMMON.MPI"
 190 !      include "COMMON.FREE"
 191       include "mpif.h"
 192       integer :: n,chunk,i,j,ii,remainder
 193       integer :: kolorW,key,ierror,errcode
 194       logical :: lprint
 195       lprint=.true.
 196 !
 197 ! Split the communicator if independent runs for different parameter
 198 ! sets will be performed.
 199 !
 200       if (nparmset.eq.1 .or. .not.separate_parset) then
 201         WHAM_COMM = MPI_COMM_WORLD
 202       else if (separate_parset) then
 203         if (nprocs.lt.nparmset) then
 204           write (iout,*) &
 205        "*** Cannot split parameter sets for fewer processors than sets",&
 206         nprocs,nparmset
 207           call MPI_Finalize(ierror)
 208           stop
 209         endif
 210         write (iout,*) "nparmset",nparmset
 211         nprocs = nprocs/nparmset
 212         kolorW = me/nprocs
 213         key = mod(me,nprocs)
 214         write (iout,*) "My old rank",me," kolor",kolorW," key",key
 215         call MPI_Comm_split(MPI_COMM_WORLD,kolorW,key,WHAM_COMM,ierror)
 216         call MPI_Comm_size(WHAM_COMM,nprocs,ierror)
 217         call MPI_Comm_rank(WHAM_COMM,me,ierror)
 218         write (iout,*) "My new rank",me," comm size",nprocs
 219         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
 220          " WHAM_COMM",WHAM_COMM
 221         myparm=kolorW+1
 222         write (iout,*) "My parameter set is",myparm
 223         call flush(iout)
 224       else
 225         myparm=nparmset
 226       endif
 227       Me1 = Me
 228       Nprocs1 = Nprocs
 229       return
 230       end subroutine proc_groups
 231 #endif
 232 !------------------------------------------------------------------------------
 233 #ifdef AIX
 234       subroutine flush(iu)
 235       call flush_(iu)
 236       return
 237       end subroutine flush
 238 #endif
 239 !-----------------------------------------------------------------------------
 240       subroutine promienie(*)
 241
 242       use io_units
 243       use names
 244       use io_base, only:ucase
 245       use energy_data, only:sigma0,dsc,dsc_inv
 246       use wham_data
 247       use w_compar_data
 248       implicit none
 249 !      include 'DIMENSIONS'
 250 !      include 'COMMON.CONTROL'
 251 !      include 'COMMON.INTERACT'
 252 !      include 'COMMON.IOUNITS'
 253 !      include 'COMMON.CONTPAR'
 254 !      include 'COMMON.LOCAL'
 255       integer ::i,j
 256       real(kind=8) :: facont=1.569D0  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
 257       character(len=8) :: contfunc
 258       character(len=8) :: contfuncid(5)=reshape((/'GB      ',&
 259                       'DIST    ','CEN     ','ODC     ','SIG     '/),shape(contfuncid))
 260 !el      character(len=8) ucase
 261       call getenv("CONTFUNC",contfunc)
 262       contfunc=ucase(contfunc)
 263       do icomparfunc=1,5
 264         if (contfunc.eq.contfuncid(icomparfunc)) goto 10
 265       enddo
 266    10 continue
 267       write (iout,*) "Sidechain contact function is ",contfunc,&
 268         "icomparfunc",icomparfunc
 269       do i=1,ntyp
 270         do j=1,ntyp
 271           if (icomparfunc.lt.3) then
 272             read(isidep1,*) chi_comp(i,j),chip_comp(i,j),sig_comp(i,j),&
 273              sc_cutoff(i,j)
 274           else if (icomparfunc.lt.5) then
 275             read(isidep1,*) sc_cutoff(i,j)
 276           else if (icomparfunc.eq.5) then
 277             sc_cutoff(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)*facont
 278           else
 279             write (iout,*) "Error - Unknown contact function"
 280             return 1
 281           endif
 282         enddo
 283       enddo
 284       close (isidep1)
 285       do i=1,ntyp1
 286         if (i.eq.10 .or. i.eq.ntyp1) then
 287           dsc_inv(i)=0.0d0
 288         else
 289           dsc_inv(i)=1.0d0/dsc(i)
 290         endif
 291       enddo
 292       return
 293       end subroutine promienie
 294 !-----------------------------------------------------------------------------
 295       subroutine alloc_wham_arrays
 296
 297       use names
 298       use geometry_data, only:nres
 299       use energy_data, only:maxcont
 300       use wham_data
 301       use w_compar_data
 302       integer :: i,j,k,l
 303 !-------------------------
 304 ! COMMON.FREE
 305 !      common /wham/
 306       allocate(stot(nslice)) !(maxslice)
 307       do i=1,nslice
 308         stot(i)=0
 309       enddo
 310       allocate(Kh(nQ,MaxR,MaxT_h,nParmSet),q0(nQ,MaxR,MaxT_h,nParmSet))!(MaxQ,MaxR,MaxT_h,max_parm)
 311       allocate(f(maxR,maxT_h,nParmSet)) !(maxR,maxT_h,max_parm)
 312       allocate(beta_h(maxT_h,nParmSet)) !(MaxT_h,max_parm)
 313       allocate(nR(maxT_h,nParmSet),nRR(maxT_h,nParmSet)) !(maxT_h,max_parm)
 314       allocate(snk(MaxR,MaxT_h,nParmSet,nSlice)) !(MaxR,MaxT_h,max_parm,MaxSlice)
 315 !      do i=1,MaxR
 316 !        do j=1,MaxT_h
 317 !          do k=1,nParmSet
 318 !            do l=1,nSlice
 319 !              snk(i,j,k,l)=0
 320 !            enddo
 321 !          enddo
 322 !        enddo
 323 !      enddo
 324
 325       allocate(totraj(maxR,nParmSet)) !(maxR,max_parm)
 326
 327       allocate(nT_h(nParmSet))!(max_parm)
 328       allocate(replica(nParmSet))
 329       allocate(umbrella(nParmSet))
 330       allocate(read_iset(nParmSet))
 331 !      allocate(nT_h(nParmSet))
 332 !-------------------------
 333 ! COMMON.PROT
 334 !      common /protein/
 335       allocate(ntot(nslice))  !(maxslice)
 336 !      allocatable :: isampl        !(max_parm)
 337 !-------------------------
 338 ! COMMON.PROTFILES
 339 !      common /protfil/
 340       allocate(protfiles(maxfile_prot,2,MaxR,MaxT_h,nParmSet)) !(maxfile_prot,2,MaxR,MaxT_h,Max_Parm)
 341       allocate(nfile_bin(MaxR,MaxT_h,nParmSet))
 342       allocate(nfile_asc(MaxR,MaxT_h,nParmSet))
 343       allocate(nfile_cx(MaxR,MaxT_h,nParmSet))
 344       allocate(rec_start(MaxR,MaxT_h,nParmSet))
 345       allocate(rec_end(MaxR,MaxT_h,nParmSet)) !(MaxR,MaxT_h,Max_Parm)
 346 !-------------------------
 347 ! COMMON.OBCINKA
 348 !      common /obcinka/
 349       allocate(time_start_collect(maxR,MaxT_h,nParmSet))
 350       allocate(time_end_collect(maxR,MaxT_h,nParmSet)) !(maxR,MaxT_h,Max_Parm)
 351 !-------------------------
 352 ! COMMON.CONTPAR
 353 !      common /contpar/
 354       allocate(sig_comp(ntyp,ntyp),chi_comp(ntyp,ntyp),&
 355         chip_comp(ntyp,ntyp),sc_cutoff(ntyp,ntyp)) !(ntyp,ntyp)
 356 !-------------------------
 357 ! COMMON.PEPTCONT
 358 !      common /peptcont/
 359       allocate(icont_pept_ref(2,maxcont)) !(2,maxcont)
 360 !      allocate(ncont_frag_ref()) !(mmaxfrag)
 361 !      allocate(icont_frag_ref(2,maxcont)) !(2,maxcont,mmaxfrag)
 362       allocate(isec_ref(nres)) !(maxres)
 363 !-------------------------
 364 ! COMMON.VAR
 365 ! Angles from experimental structure
 366 !      common /varref/
 367       allocate(vbld_ref(nres),theta_ref(nres),&
 368         phi_ref(nres),alph_ref(nres),omeg_ref(nres)) !(maxres)
 369 !-------------------------
 370       end subroutine alloc_wham_arrays
 371 !-----------------------------------------------------------------------------
 372 !-----------------------------------------------------------------------------