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,hpb_partition
  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       endif
 133       if (constr_dist.gt.0) then
 134          call read_dist_constr
 135          call hpb_partition
 136       endif
 137       write (iout,*) "Begin read_database"
 138       call flush(iout)
 139       call read_database(*10)
 140       write (iout,*) "Finished read_database"
 141       call flush(iout)
 142       if (separate_parset) nparmset=1
 143       do islice=1,nslice
 144         if (ntot(islice).gt.0) then
 145 #ifdef MPI
 146           call work_partition(islice,.true.)
 147           write (iout,*) "work_partition OK"
 148           call flush(iout)
 149 #endif
 150           write (iout,*) "call enecalc",islice,nslice
 151           call enecalc(islice,*10)
 152           write (iout,*) "enecalc OK"
 153           call flush(iout)
 154           call WHAMCALC(islice,*10)
 155           write (iout,*) "wham_calc OK"
 156           call flush(iout)
 157           call write_dbase(islice,*10)
 158           write (iout,*) "write_dbase OK"
 159           call flush(iout)
 160           if (ensembles.gt.0) then
 161             call make_ensembles(islice,*10)
 162             write (iout,*) "make_ensembles OK"
 163             call flush(iout)
 164           endif
 165         endif
 166       enddo
 167 #ifdef MPI
 168       call MPI_Finalize( IERROR )
 169 #endif
 170       stop
 171    10 write (iout,*) "Error termination of the program"
 172 #ifdef MPI
 173       call MPI_Finalize( IERROR )
 174 #endif
 175       stop
 176       end program wham_multparm
 177 !------------------------------------------------------------------------------
 178 !
 179 !------------------------------------------------------------------------------
 180 #ifdef MPI
 181       subroutine proc_groups
 182 ! Split the processors into the Master and Workers group, if needed.
 183       use io_units
 184       use MPI_data
 185       use wham_data
 186
 187       implicit none
 188 !      include "DIMENSIONS"
 189 !      include "DIMENSIONS.ZSCOPT"
 190 !      include "DIMENSIONS.FREE"
 191 !      include "COMMON.IOUNITS"
 192 !      include "COMMON.MPI"
 193 !      include "COMMON.FREE"
 194       include "mpif.h"
 195       integer :: n,chunk,i,j,ii,remainder
 196       integer :: kolorW,key,ierror,errcode
 197       logical :: lprint
 198       lprint=.true.
 199 !
 200 ! Split the communicator if independent runs for different parameter
 201 ! sets will be performed.
 202 !
 203       if (nparmset.eq.1 .or. .not.separate_parset) then
 204         WHAM_COMM = MPI_COMM_WORLD
 205       else if (separate_parset) then
 206         if (nprocs.lt.nparmset) then
 207           write (iout,*) &
 208        "*** Cannot split parameter sets for fewer processors than sets",&
 209         nprocs,nparmset
 210           call MPI_Finalize(ierror)
 211           stop
 212         endif
 213         write (iout,*) "nparmset",nparmset
 214         nprocs = nprocs/nparmset
 215         kolorW = me/nprocs
 216         key = mod(me,nprocs)
 217         write (iout,*) "My old rank",me," kolor",kolorW," key",key
 218         call MPI_Comm_split(MPI_COMM_WORLD,kolorW,key,WHAM_COMM,ierror)
 219         call MPI_Comm_size(WHAM_COMM,nprocs,ierror)
 220         call MPI_Comm_rank(WHAM_COMM,me,ierror)
 221         write (iout,*) "My new rank",me," comm size",nprocs
 222         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
 223          " WHAM_COMM",WHAM_COMM
 224         myparm=kolorW+1
 225         write (iout,*) "My parameter set is",myparm
 226         call flush(iout)
 227       else
 228         myparm=nparmset
 229       endif
 230       Me1 = Me
 231       Nprocs1 = Nprocs
 232       return
 233       end subroutine proc_groups
 234 #endif
 235 !------------------------------------------------------------------------------
 236 #ifdef AIX
 237       subroutine flush(iu)
 238       call flush_(iu)
 239       return
 240       end subroutine flush
 241 #endif
 242 !-----------------------------------------------------------------------------
 243       subroutine promienie(*)
 244
 245       use io_units
 246       use names
 247       use io_base, only:ucase
 248       use energy_data, only:sigma0,dsc,dsc_inv
 249       use wham_data
 250       use w_compar_data
 251       implicit none
 252 !      include 'DIMENSIONS'
 253 !      include 'COMMON.CONTROL'
 254 !      include 'COMMON.INTERACT'
 255 !      include 'COMMON.IOUNITS'
 256 !      include 'COMMON.CONTPAR'
 257 !      include 'COMMON.LOCAL'
 258       integer ::i,j
 259       real(kind=8) :: facont=1.569D0  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
 260       character(len=8) :: contfunc
 261       character(len=8) :: contfuncid(5)=reshape((/'GB      ',&
 262                       'DIST    ','CEN     ','ODC     ','SIG     '/),shape(contfuncid))
 263 !el      character(len=8) ucase
 264       call getenv("CONTFUNC",contfunc)
 265       contfunc=ucase(contfunc)
 266       do icomparfunc=1,5
 267         if (contfunc.eq.contfuncid(icomparfunc)) goto 10
 268       enddo
 269    10 continue
 270       write (iout,*) "Sidechain contact function is ",contfunc,&
 271         "icomparfunc",icomparfunc
 272       do i=1,ntyp
 273         do j=1,ntyp
 274           if (icomparfunc.lt.3) then
 275             read(isidep1,*) chi_comp(i,j),chip_comp(i,j),sig_comp(i,j),&
 276              sc_cutoff(i,j)
 277           else if (icomparfunc.lt.5) then
 278             read(isidep1,*) sc_cutoff(i,j)
 279           else if (icomparfunc.eq.5) then
 280             sc_cutoff(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)*facont
 281           else
 282             write (iout,*) "Error - Unknown contact function"
 283             return 1
 284           endif
 285         enddo
 286       enddo
 287       close (isidep1)
 288       do i=1,ntyp1
 289         if (i.eq.10 .or. i.eq.ntyp1) then
 290           dsc_inv(i)=0.0d0
 291         else
 292           dsc_inv(i)=1.0d0/dsc(i)
 293         endif
 294       enddo
 295       return
 296       end subroutine promienie
 297 !-----------------------------------------------------------------------------
 298       subroutine alloc_wham_arrays
 299
 300       use names
 301       use geometry_data, only:nres
 302       use energy_data, only:maxcont
 303       use wham_data
 304       use w_compar_data
 305       integer :: i,j,k,l
 306 !-------------------------
 307 ! COMMON.FREE
 308 !      common /wham/
 309       allocate(stot(nslice)) !(maxslice)
 310       do i=1,nslice
 311         stot(i)=0
 312       enddo
 313       allocate(Kh(nQ,MaxR,MaxT_h,nParmSet),q0(nQ,MaxR,MaxT_h,nParmSet))!(MaxQ,MaxR,MaxT_h,max_parm)
 314       allocate(f(maxR,maxT_h,nParmSet)) !(maxR,maxT_h,max_parm)
 315       allocate(beta_h(maxT_h,nParmSet)) !(MaxT_h,max_parm)
 316       allocate(nR(maxT_h,nParmSet),nRR(maxT_h,nParmSet)) !(maxT_h,max_parm)
 317       allocate(snk(MaxR,MaxT_h,nParmSet,nSlice)) !(MaxR,MaxT_h,max_parm,MaxSlice)
 318 !      do i=1,MaxR
 319 !        do j=1,MaxT_h
 320 !          do k=1,nParmSet
 321 !            do l=1,nSlice
 322 !              snk(i,j,k,l)=0
 323 !            enddo
 324 !          enddo
 325 !        enddo
 326 !      enddo
 327
 328       allocate(totraj(maxR,nParmSet)) !(maxR,max_parm)
 329
 330       allocate(nT_h(nParmSet))!(max_parm)
 331       allocate(replica(nParmSet))
 332       allocate(umbrella(nParmSet))
 333       allocate(read_iset(nParmSet))
 334 !      allocate(nT_h(nParmSet))
 335 !-------------------------
 336 ! COMMON.PROT
 337 !      common /protein/
 338       allocate(ntot(nslice))  !(maxslice)
 339 !      allocatable :: isampl        !(max_parm)
 340 !-------------------------
 341 ! COMMON.PROTFILES
 342 !      common /protfil/
 343       allocate(protfiles(maxfile_prot,2,MaxR,MaxT_h,nParmSet)) !(maxfile_prot,2,MaxR,MaxT_h,Max_Parm)
 344       allocate(nfile_bin(MaxR,MaxT_h,nParmSet))
 345       allocate(nfile_asc(MaxR,MaxT_h,nParmSet))
 346       allocate(nfile_cx(MaxR,MaxT_h,nParmSet))
 347       allocate(rec_start(MaxR,MaxT_h,nParmSet))
 348       allocate(rec_end(MaxR,MaxT_h,nParmSet)) !(MaxR,MaxT_h,Max_Parm)
 349 !-------------------------
 350 ! COMMON.OBCINKA
 351 !      common /obcinka/
 352       allocate(time_start_collect(maxR,MaxT_h,nParmSet))
 353       allocate(time_end_collect(maxR,MaxT_h,nParmSet)) !(maxR,MaxT_h,Max_Parm)
 354 !-------------------------
 355 ! COMMON.CONTPAR
 356 !      common /contpar/
 357       allocate(sig_comp(ntyp,ntyp),chi_comp(ntyp,ntyp),&
 358         chip_comp(ntyp,ntyp),sc_cutoff(ntyp,ntyp)) !(ntyp,ntyp)
 359 !-------------------------
 360 ! COMMON.PEPTCONT
 361 !      common /peptcont/
 362       allocate(icont_pept_ref(2,maxcont)) !(2,maxcont)
 363 !      allocate(ncont_frag_ref()) !(mmaxfrag)
 364 !      allocate(icont_frag_ref(2,maxcont)) !(2,maxcont,mmaxfrag)
 365       allocate(isec_ref(nres)) !(maxres)
 366 !-------------------------
 367 ! COMMON.VAR
 368 ! Angles from experimental structure
 369 !      common /varref/
 370       allocate(vbld_ref(nres),theta_ref(nres),&
 371         phi_ref(nres),alph_ref(nres),omeg_ref(nres)) !(maxres)
 372 !-------------------------
 373       end subroutine alloc_wham_arrays
 374 !-----------------------------------------------------------------------------
 375 !-----------------------------------------------------------------------------