source/maxlik/src-Fmatch_safe/make_list_MD.F

   1       subroutine make_list_MD(lprn)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       integer IERROR,ErrCode,Status(MPI_STATUS_SIZE,10)
   8       integer req(10),msg_in(5),msg_out(5),address,bufsize
   9       character*1 buffer(8*(maxstr_proc*maxres6*max_ene+8000))
  10       include "COMMON.MPI"
  11 #endif
  12       include "COMMON.WEIGHTS"
  13       include "COMMON.WEIGHTDER"
  14       include "COMMON.ENERGIES"
  15       include "COMMON.IOUNITS"
  16       include "COMMON.VMCPAR"
  17       include "COMMON.NAMES"
  18       include "COMMON.INTERACT"
  19       include "COMMON.TIME1"
  20       include "COMMON.CHAIN"
  21       include "COMMON.PROTFILES"
  22       include "COMMON.VAR"
  23       include "COMMON.GEO"
  24       include "COMMON.OPTIM"
  25       include "COMMON.CLASSES"
  26       include "COMMON.COMPAR"
  27       include "COMMON.TORSION"
  28       include "COMMON.FMATCH"
  29       include "COMMON.DERIV"
  30       include "COMMON.PROTNAME"
  31 C Define local variables
  32       integer i,ii,iii,ibatch,kkk,jj,j,k,kk,l,iprot,ib,nn
  33       integer ipass_conf,istart_conf,iend_conf,Previous,Next
  34       double precision energia(0:max_ene)
  35       double precision tcpu_ini,tcpu_fin,tcpu
  36       integer iroof,icant
  37       external iroof,icant
  38       logical lprn,first_call
  39       logical*1 lflag(maxstr)
  40       integer list_conf_MD_(maxstr,maxprot)
  41       double precision ftune_eps
  42       integer ilen
  43       external ilen
  44 #ifdef MPI
  45 c      double precision e_total_(maxstr_proc)
  46 #endif
  47       write (iout,*) "Making the worklist of MD conformations."
  48       call flush(iout)
  49       calc_grad=.true.
  50       tcpu_ini = tcpu()
  51       do iprot=1,nprot
  52         do i=1,nconf_forc(iprot)
  53           list_conf_MD(i,iprot)=i
  54         enddo
  55         ntot_work_MD(iprot)=nconf_forc(iprot)
  56         write (iout,*) "iprot",iprot," ntot_work_MD",ntot_work_MD(iprot)
  57         call flush(iout)
  58       enddo
  59 #ifdef MPI
  60 c
  61 c Divide the whole database between processors
  62 c
  63       Previous = me1-1
  64       Next = me1+1
  65       if (Previous.lt.0) Previous = Nprocs-1
  66       if (Next.ge.Nprocs) Next = 0
  67       write (iout,*) "Make_list_MD: first call of work_partition"
  68       call flush(iout)
  69       call work_partition_MD(lprn)
  70       write (iout,*) "Make_list_MD: end first call of work_partition"
  71       call flush(iout)
  72 #endif
  73 c
  74 c Loop over proteins
  75 c
  76       DO iprot=1,nprot
  77
  78        if (.not.fmatch(iprot)) cycle
  79        call restore_molinfo(iprot)
  80 c Loop over the conformations of protein IPROT assigned to the current processor.
  81 c The conformations are read off a DA scratchfile and processed in passes, each
  82 c of which requires no more than MAXSTR_PROC conformations to be stored in memory
  83 c simultaneously.
  84 c
  85       write (iout,*) "iprot",iprot
  86       write (iout,*) "Opening ",bprotfiles_MD(iprot)
  87       call flush(iout)
  88       open (icbase,file=bprotfiles_MD(iprot),status="old",
  89      &    form="unformatted",access="direct",recl=lenrec_MD(iprot))
  90       write (iout,*) "Opening ",bforcefiles(iprot)
  91       call flush(iout)
  92       if (.not.mod_other_params)
  93      &  open (ientout,file=bforcefiles(iprot),status="unknown",
  94      &    form="unformatted",access="direct",recl=lenrec_forces(iprot))
  95       write (iout,*) bforcefiles(iprot)," opened"
  96
  97 #ifdef MPI
  98       nchunk_conf_MD(iprot) = iroof(scount_MD(me1,iprot),maxstr_proc)
  99       write (iout,*)"Protein",iprot," from MD energy evaluation in",
 100      &   nchunk_conf_MD(iprot)," passes."
 101       call flush(iout)
 102       ipass_conf=0
 103       jj=0
 104       nn=0
 105       do i=indstart_MD(me1,iprot),indend_MD(me1,iprot),maxstr_proc
 106       ipass_conf=ipass_conf+1
 107       write (iout,*) "MAKE_LIST_MD: Pass",ipass_conf
 108       istart_conf=i
 109       iend_conf=min0(i+maxstr_proc-1,indend_MD(me1,iprot))
 110       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
 111 #else
 112       nchunk_ene_MD(iprot) = iroof(ntot_MD(iprot),maxstr)
 113       nchunk_conf_MD(iprot) = iroof(ntot_MD(iprot),maxstr_proc)
 114       write (iout,*)"Protein",iprot," energy evaluation in",
 115      &   nchunk_conf_MD," passes."
 116       ipass_conf=0
 117
 118       do i=1,ntot_MD(iprot),maxstr_proc
 119
 120       ipass_conf=ipass_conf+1
 121       write (iout,*) "MAKE_LIST: Pass",ipass_conf_MD
 122       istart_conf=i
 123       iend_conf=min0(i+maxstr_proc-1,ntot_MD(iprot))
 124 #endif
 125       ii=0
 126 c Read the chunk of conformations off a DA scratchfile.
 127 c
 128       call daread_MDcoords(iprot,istart_conf,iend_conf)
 129 c
 130 c Compute the energies of the conformations currently in memory and compute
 131 c the lowest energies.
 132 c
 133       do iii=istart_conf,iend_conf
 134         ii=ii+1
 135         call restore_MDcoords(iprot,ii)
 136         call int_from_cart1(.false.)
 137 #ifdef DEBUG
 138         write (iout,*) "Makelist_MD: Before etotal",iii,i
 139         call flush(iout)
 140 #endif
 141         call zerograd
 142         call etotal(energia(0))
 143         call cartgrad
 144 #ifdef DEBUG
 145         write (iout,*) "After etotal",iii
 146         call flush(iout)
 147         call enerprint(energia(0))
 148 c          write (iout,'(i5,16(1pe12.4))') i,
 149 c     &    (energia(j),j=1,n_ene),energia(0)
 150         call flush(iout)
 151 #endif
 152 #ifdef DEBUG
 153         write (iout,*) "Conformation:",iii
 154         write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 155         write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 156         write (iout,'(8f10.4)') (vbld(k),k=2,nres)
 157         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
 158         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 159         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 160         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 161         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 162         call enerprint(energia(0))
 163 #endif
 164 c        e_total_MD(iii,iprot)=energia(0)
 165         do k=1,nres
 166           do l=1,3
 167             do j=1,n_ene
 168               forctb(j,l,k,ii,iprot)=-gcompon(j,l,k)
 169             enddo
 170           enddo
 171         enddo
 172         kk=nres
 173         do k=nnt,nct
 174           if (itype(k).eq.10 .or. itype(k).eq.ntyp1) cycle
 175           kk=kk+1
 176           do l=1,3
 177             do j=1,n_ene
 178               forctb(j,l,kk,ii,iprot)=-gcomponx(j,l,k)
 179             enddo
 180           enddo
 181         enddo
 182 #ifdef DEBUG
 183         write (iout,'(i5,20(1pe12.4))') iii,
 184      &    (energia(j),j=1,n_ene),energia(0)
 185         call flush(iout)
 186 #endif
 187         if (energia(0).ge.1.0d7) then
 188           write (iout,*)
 189      &  "MAKE_LIST_MD: energy too high for protein",iprot,
 190      &  " point",iii,ii,", conformationn skipped."
 191           write (iout,*) "The components of the energy are:"
 192           call enerprint(energia(0))
 193           write (iout,*) "Conformation:",iii,ii
 194           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 195           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 196           write (iout,'(8f10.4)') (vbld(k),k=2,nres)
 197           write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
 198           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 199           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 200           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 201           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 202         else
 203           nn=nn+1
 204           list_conf_MD_(nn,iprot)=iii
 205         endif
 206       enddo ! iii
 207 C Change AL 12/30/2017
 208 c
 209       write (iout,*) "mod_other_params ",mod_other_params
 210       if (.not. mod_other_params) then
 211         write (iout,*) "make_list_MD: calling dawrite_forces"
 212         write (iout,*) "istart_conf",istart_conf,
 213      &    " iend_conf",iend_conf
 214
 215         call dawrite_forces(iprot,istart_conf,iend_conf,ientout)
 216 c        write (iout,*) "Finish dawrite_forces"
 217 c        call flush(iout)
 218 #ifdef MPI
 219 c Distribute force components through ring
 220 c        write (iout,*) "Calling MPI_Barrier"
 221 c        call flush(iout)
 222         call MPI_Barrier(WHAM_COMM,IERROR)
 223         write (iout,*) "Processes synchronized in make_list_MD"
 224         call flush(iout)
 225         msg_out(1)= 5*me1+1000*ipass_conf+6
 226         msg_out(2)= 5*me1+1000*ipass_conf+7
 227         msg_out(3)= 5*me1+1000*ipass_conf+8
 228         msg_out(4)= 5*me1+1000*ipass_conf+9
 229         msg_out(5)= 5*me1+1000*ipass_conf+10
 230         do iii=1,Nprocs-1
 231 c Send the current force tables to the right neighbor
 232 c Receive the force tables produced by processor kkk from the left neighbor
 233           kkk = mod(me1-iii+NProcs,Nprocs)
 234           msg_in(1)= 5*kkk+1000*ipass_conf+6
 235           msg_in(2)= 5*kkk+1000*ipass_conf+7
 236           msg_in(3)= 5*kkk+1000*ipass_conf+8
 237           msg_in(4)= 5*kkk+1000*ipass_conf+9
 238           msg_in(5)= 5*kkk+1000*ipass_conf+10
 239           write (iout,*) "me1",me1," iii",iii," Previous",Previous,
 240      &     " Next",Next," kkk",kkk
 241           write (iout,*) "msg_in",msg_in
 242           write (iout,*) "msg_out",msg_out
 243           call flush(iout)
 244           write (iout,*) "Processor",me1," Start Send and receive"
 245           call flush(iout)
 246           call MPI_Send(istart_conf,1,MPI_INTEGER,Next,msg_out(1),
 247      &     WHAM_COMM,IERROR)
 248           write (iout,*) "Send",msg_out(1)," complete"
 249           call flush(iout)
 250           call MPI_Recv(istart_conf,1,MPI_INTEGER,Previous,
 251      &     msg_in(1),WHAM_COMM,STATUS(1,6),IERROR)
 252           write (iout,*) "Recv",msg_in(1)," complete"
 253           call flush(iout)
 254           call MPI_Send(iend_conf,1,MPI_INTEGER,Next,msg_out(2),
 255      &     WHAM_COMM,IERROR)
 256           write (iout,*) "Send",msg_out(2)," complete"
 257           call flush(iout)
 258           call MPI_Recv(iend_conf,1,MPI_INTEGER,Previous,
 259      &     msg_in(2),WHAM_COMM,STATUS(1,7),IERROR)
 260           write (iout,*) "Recv",msg_in(2)," complete"
 261           call flush(iout)
 262           call MPI_Buffer_Attach(buffer(1),
 263      &     8*(scount_MD(Me1,iprot)*maxres6*max_ene+800),
 264      &     IERROR)
 265           call MPI_BSend(forctb(1,1,1,1,iprot),
 266      &          scount_MD(me1,iprot)*maxres6*max_ene,
 267      &          MPI_DOUBLE_PRECISION,Next,msg_out(3),
 268      &          WHAM_COMM,IERROR)
 269           write (iout,*) "Send",msg_out(3)," complete (forctb)"
 270           call flush(iout)
 271           call MPI_Recv(forctb(1,1,1,1,iprot),
 272      &          scount_MD(Previous,iprot)*maxres6*max_ene,
 273      &          MPI_DOUBLE_PRECISION,Previous,msg_in(3),WHAM_COMM,
 274      &          STATUS(1,8),IERROR)
 275           write (iout,*) "Recv",msg_in(3)," complete (forctb)"
 276           call flush(iout)
 277           call MPI_Buffer_Detach(address,bufsize,IERROR)
 278           write (iout,*) "Comm with proc",iii," istart_conf",
 279      &      istart_conf," iend_conf",iend_conf
 280           write (iout,*) "Calling dawrite_forces"
 281           call dawrite_forces(iprot,istart_conf,iend_conf,ientout)
 282           do k=1,5
 283             msg_out(k)=msg_in(k)
 284           enddo
 285         enddo
 286 #endif
 287         if (.not.mod_other_params) close(ientout)
 288       endif
 289
 290       enddo
 291       close (icbase)
 292       ntot_work_MD(iprot)=nn
 293
 294       ENDDO ! iprot
 295
 296 c Divide the work again based on the current lists of conformations
 297       DO IPROT=1,NPROT
 298 #ifdef DEBUG
 299        write (iout,*) "Processor",me," loop IPROT",iprot
 300        write (iout,*) "fmatch ",fmatch(iprot)
 301 #endif
 302        if (.not.fmatch(iprot)) cycle
 303 #ifdef MPI
 304 c
 305 c All workers get the complete list of conformations.
 306 c
 307         call MPI_Allgather(ntot_work_MD(iprot),1,MPI_INTEGER,
 308      &    scount_MD(0,iprot),1,MPI_INTEGER,Comm1,IERROR)
 309         idispl_MD(0,iprot)=0
 310         do i=1,nprocs1-1
 311           idispl_MD(i,iprot)=idispl_MD(i-1,iprot)+scount_MD(i-1,iprot)
 312         enddo
 313 #ifdef DEBUG
 314         write (iout,*) "Protein",iprot," MD Scount and Idispl"
 315         do i=0,nprocs1-1
 316           write (iout,*) i,scount_MD(i,iprot),idispl_MD(i,iprot)
 317         enddo
 318         write (iout,*) "Protein",iprot,
 319      &    " local list of MD conformations of processor",me
 320         do i=1,ntot_work_MD(iprot)
 321           write(iout,*) i,list_conf_MD(i,iprot)
 322         enddo
 323         write (iout,*) "Before REDUCE: ntot_work_MD",ntot_work_MD(iprot)
 324         call flush(iout)
 325 #endif
 326         call MPI_Allreduce(ntot_work_MD(iprot),nn,1,MPI_INTEGER,MPI_SUM,
 327      &    Comm1,IERROR)
 328         ntot_work_MD(iprot)=nn
 329 #ifdef DEBUG
 330         write (iout,*) "After REDUCE: ntot_work_MD",ntot_work_MD(iprot)
 331         call flush(iout)
 332 #endif
 333         call MPI_Allgatherv(list_conf_MD_(1,iprot),
 334      &    scount_MD(me1,iprot),MPI_INTEGER,list_conf_MD(1,iprot),
 335      &   scount_MD(0,iprot),idispl_MD(0,iprot),MPI_INTEGER,Comm1,IERROR)
 336 #ifdef DEBUG
 337         write (iout,*) "Protein",iprot,
 338      &    " global list of MD conformations of processor",me
 339         do i=1,ntot_work_MD(iprot)
 340           write(iout,'(2i5)')i,list_conf_MD(i,iprot)
 341 c     &      ,(nu(k,i,iprot),k=1,natlike(iprot))
 342         enddo
 343         call flush(iout)
 344 #endif
 345 #else
 346         do i=1,ntot_work_MD(iprot)
 347           list_conf_MD(i,iprot)=list_conf_MD_(i,iprot)
 348         enddo
 349 #endif
 350       ENDDO !  IPROT
 351 c
 352       call work_partition_MD(.true.)
 353 c
 354       do iprot=1,nprot
 355        if (.not.fmatch(iprot)) cycle
 356 #ifdef MPI
 357         do i=1,ntot_work_MD(iprot)
 358           i2ii_MD(i,iprot)=0
 359         enddo
 360         ii=0
 361         do i=indstart_MD(me1,iprot),indend_MD(me1,iprot)
 362           ii=ii+1
 363           i2ii_MD(i,iprot)=ii
 364         enddo
 365 #else
 366         do i=1,ntot_work_MD(iprot)
 367           i2ii_MD(i,iprot)=i
 368         endif
 369 #endif
 370 c
 371         tsil_MD=0
 372         do i=1,ntot_work_MD(iprot)
 373           tsil_MD(list_conf_MD(i,iprot),iprot)=i
 374         enddo
 375 #ifdef DEBUG
 376         write (iout,*) "Protein",i," List-to-conformation mapping"
 377         do i=1,ntot_work_MD(iprot)
 378           write(iout,*) i,tsil_MD(i,iprot)
 379         enddo
 380 #endif
 381       enddo       ! iprot
 382
 383       write (iout,*)
 384       write (iout,*)
 385      &  "Numbers of MD conformations after applying the cutoff"
 386       write (iout,*)
 387       do iprot=1,nprot
 388         if (.not.fmatch(iprot)) cycle
 389         write (iout,'(a,2x,a,i10)') "Protein",
 390      &    protname(iprot)(:ilen(protname(iprot))),
 391      &    ntot_work_MD(iprot)
 392       enddo
 393 c
 394
 395       tcpu_fin = tcpu() - tcpu_ini
 396       write (iout,*) "Time for creating list of MD conformations",
 397      &  tcpu_fin
 398       call flush(iout)
 399       t_func = t_func + tcpu_fin
 400       return
 401       end