source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/pmfread.F

   1       subroutine pmfread(*)
   2       include "DIMENSIONS"
   3       include "DIMENSIONS.ZSCOPT"
   4 #ifdef MPI
   5       include 'mpif.h'
   6       include 'COMMON.MPI'
   7       include 'COMMON.MPISET'
   8       integer ierror
   9       double precision bufin(maxpoint),bufout(maxpoint)
  10 #endif
  11       include "COMMON.PMFDERIV"
  12       include "COMMON.IOUNITS"
  13       include "COMMON.WEIGHTS"
  14 c      include "COMMON.MARQ"
  15 c      character*64 fname(8)
  16 c      common /kuku/ it,jt
  17       character*1 typ(-2:2) /'p','a','G','A','P'/
  18       character*64 torfile
  19       character*16 pmf_print
  20       integer ilen
  21       external ilen
  22       logical lprn
  23       double precision rjunk,sumw(4)
  24       write (iout,*) "Calling PMFREAD",torsion_pmf,turn_pmf,eello_pmf
  25       call getenv("PMF_DIR",pmf_dir)
  26       call getenv("PMF_INFIX",pmf_infix)
  27       write (iout,*) "PMF_DIR ",PMF_DIR
  28       write (iout,*) "PMF_INFIX",PMF_INFIX
  29       call getenv("PMF_PRINT",pmf_print)
  30       lprn=index(pmf_print,"PMF_PRINT").gt.0
  31 c      logical lder
  32 c      do i=1,n
  33 c        read(iparin,'(e20.10,1x,a)') x(i),parname(i)
  34 c        print *,i,x(i)," ",parname(i)
  35 c      enddo
  36 c      maxit=100
  37 c      maxmar=10
  38 c      rl=1.0d2
  39 c      vmarq=1.0d1
  40 c      tolx=1.0d-3
  41 c      tolf=1.0d-2
  42 c      rtolf=1.0d-3
  43 c      tollam=1.0d2
  44 c      read (iparin,*,end=12,err=12)
  45 c     & maxit,maxmar,rl,vmarq,tolx,tolf,rtolf,tollam
  46 c   12 continue
  47 c      write (iout,'(//21(1h*)/a/21(1h*))') "Initial parameters:"
  48 c      do i=1,n
  49 c        write (iout,'(i3,1x,a16,f10.5)') i,parname(i),x(i)
  50 c      enddo
  51
  52       ifunmode=1
  53       m=3
  54 c      write (iout,'(a,i3)') 'number of variables  ',m
  55 c      write (iout,'(a,i3)') 'number of parameters ',n
  56       ii=1
  57       i1=1
  58
  59 c      close(inp)
  60
  61       npoint=0
  62
  63       do it1=0,2
  64       do it2=-2,2
  65
  66 c      do it1=1,1
  67 c      do it2=1,1
  68
  69 c      goto 1111
  70
  71       if (torsion_pmf .or. turn_pmf) then
  72
  73 #ifdef MPI
  74       write (iout,*) "Master",Master
  75       if (me.eq.Master) then
  76 #endif
  77       torfile = PMF_DIR(:ilen(PMF_DIR))//'/'//typ(it1)//typ(it2)//"_"
  78      &  //PMF_INFIX(:ilen(PMF_INFIX))
  79      & //".eturn3"
  80       write (iout,*) "Reading torsional and turn PMFs"
  81       write (iout,'(2a)') "Torsional and turn3 potential file ",torfile
  82       open(89,file=torfile,status="old",err=99)
  83       iii = 1
  84       i1 = ii
  85       do i=1,maxpoint
  86         read (89,*,end=13,err=13) zz(1,ii),zz(2,ii),zz(3,ii),y(1,ii),
  87      &    rjunk,y(2,ii)
  88         if (y(1,ii).ge.999.0 .or. y(2,ii).ge.999.0) cycle
  89 c        if (zz(1,ii).ne.75.0d0 .or. zz(2,ii).ne.130.0d0 .or.
  90 c     &    zz(3,ii).ne.-45.0d0) cycle
  91         if (lprn)
  92      &    write(iout,'(i7,3f7.1,2f10.5)')ii,zz(1,ii),zz(2,ii),zz(3,ii),
  93      &     y(1,ii),y(2,ii)
  94         ii=ii+1
  95         iii = iii+1
  96       enddo
  97    13 np(it1,it2)=iii-1
  98       npoint=npoint+2*np(it1,it2)
  99       close (89)
 100       write (iout,*)np(it1,it2)," data points found in file"
 101       write (iout,*)npoint," data points found so far"
 102       call flush(iout)
 103       sumw=0.0d0
 104       ymin=y(1,i1)
 105       ymax=y(1,i1)
 106       do i=i1,ii-1
 107         if (y(1,i).gt.ymax) ymax=y(1,i)
 108         if (y(1,i).lt.ymin) ymin=y(1,i)
 109       enddo
 110       do i=i1,ii-1
 111         if (ifunmode.eq.0) then
 112           w(1,i)=dexp(-0.169*(y(1,i)-ymin))
 113           w(2,i)=1.0*dexp(-0.168*y(2,i))
 114         else
 115           w(1,i)=1.0d0
 116           w(2,i)=1.0d0
 117           if (y(1,i).ge.999.0) w(1,i)=0.0d0
 118           if (y(2,i).ge.999.0) w(2,i)=0.0d0
 119         endif
 120         do j=1,2
 121           sumw(j)=sumw(j)+w(j,i)
 122         enddo
 123       enddo
 124       do i=i1,ii-1
 125         do j=1,2
 126           w(j,i)=w(j,i)/sumw(j)
 127         enddo
 128       enddo
 129       close(89)
 130 #ifdef MPI
 131       np_=np(it1,it2)
 132       call work_partition_pmf(np_,.true.)
 133
 134       endif
 135
 136       write (iout,*) "After scatter"
 137       write (iout,*) "scount_pmf",(scount_pmf(i),i=1,procs)
 138       call flush(iout)
 139       call MPI_Scatter(scount_pmf,1,MPI_INTEGER,np(it1,it2),1,
 140      &     MPI_INTEGER,Master,WHAM_COMM,ierror)
 141       write (iout,*) "After scatter scount np",np(it1,it2)
 142       call flush(iout)
 143       do j=1,2
 144         do i=i1,ii-1
 145           bufin(i)=y(j,i)
 146         enddo
 147         call MPI_Scatterv(bufin(i1),scount_pmf,idispl_pmf,
 148      &    MPI_DOUBLE_PRECISION,
 149      &    bufout(1),np(it1,it2),MPI_DOUBLE_PRECISION,Master,WHAM_COMM,
 150      &    ierror)
 151         do i=1,np(it1,it2)
 152           y(j,i+i1-1)=bufout(i)
 153         enddo
 154         write (iout,*) "After scatter y",j
 155         call flush(iout)
 156         do i=i1,ii-1
 157           bufin(i)=w(j,i)
 158         enddo
 159         call MPI_Scatterv(bufin(i1),scount_pmf,idispl_pmf,
 160      &    MPI_DOUBLE_PRECISION,
 161      &    bufout(1),np(it1,it2),MPI_DOUBLE_PRECISION,Master,WHAM_COMM,
 162      &    ierror)
 163         do i=1,np(it1,it2)
 164           w(j,i+i1-1)=bufout(i)
 165         enddo
 166         write (iout,*) "After scatter w",j
 167         call flush(iout)
 168       enddo
 169       do j=1,4
 170         do i=i1,ii-1
 171           bufin(i)=zz(j,i)
 172         enddo
 173         call MPI_Scatterv(bufin(i1),scount_pmf,idispl_pmf,
 174      &    MPI_DOUBLE_PRECISION,
 175      &    bufout(1),np(it1,it2),MPI_DOUBLE_PRECISION,Master,WHAM_COMM,
 176      &    ierror)
 177         write (iout,*) "After scatter z",j
 178         call flush(iout)
 179         do i=1,np(it1,it2)
 180           zz(j,i+i1-1)=bufout(i)
 181         enddo
 182       enddo
 183       if (lprn) then
 184       write (iout,*) "Distributed etor and eturn"
 185       do i=i1,i1+np(it1,it2)-1
 186         write(iout,'(i7,3f7.1,2f10.5)')i,(zz(j,i),j=1,3),(y(j,i),j=1,2)
 187       enddo
 188       endif
 189       ii=i1+np(it1,it2)
 190       i1=ii
 191 #endif
 192 1111  continue
 193       write (iout,*)np(it1,it2)," data points found in file"
 194
 195       endif ! torsion_pmf .or. turn_pmf
 196
 197       if (eello_pmf) then
 198
 199 #define PIZDUN
 200 #ifdef PIZDUN
 201 #ifdef MPI
 202       if (me.eq.Master) then
 203 #endif
 204       torfile = PMF_DIR(:ilen(PMF_DIR))//'/'//typ(it1)//typ(it2)//"_"//
 205      &  PMF_INFIX(:ilen(PMF_INFIX))
 206      & //".eello3"
 207       write (iout,*) "Reading eello3 PMFs"
 208       write (iout,'(2a)') "local-electrostatic PMF file ",torfile
 209       open(89,file=torfile,status="old",err=99)
 210       iii = 1
 211       i1 = ii
 212       do i=1,maxpoint
 213         read (89,*,end=133,err=133) (zz(j,ii),j=1,8),(y(j,ii),j=1,4)
 214         if (y(1,ii).ge.999.0 .or. y(2,ii).ge.999.0 .or. y(3,ii).ge.999.0
 215      &     .or. y(4,ii).ge.999.0) cycle
 216         if (lprn) write(iout,'(i7,8f7.1,4f10.5)')ii,(zz(j,ii),j=1,8),
 217      &     (y(j,ii),j=1,4)
 218         ii=ii+1
 219         iii = iii+1
 220       enddo
 221   133 npel(it1,it2)=iii-1
 222       close (89)
 223       npoint=npoint+npel(it1,it2)*4
 224       sumw=0.0d0
 225       do i=i1,ii
 226         do j=1,4
 227           if (ifunmode.eq.0) then
 228             w(j,ii)=10*dexp(-0.169*(y(j,ii)))
 229           else
 230             w(j,ii)=1.0d0
 231           endif
 232         enddo
 233         do j=1,4
 234           sumw(j)=sumw(j)+w(j,i)
 235         enddo
 236       enddo
 237       do i=i1,ii-1
 238         do j=1,4
 239           w(j,i)=w(j,i)/sumw(j)
 240         enddo
 241       enddo
 242 #ifdef MPI
 243       np_=npel(it1,it2)
 244       call work_partition_pmf(np_,.true.)
 245
 246       endif
 247
 248       write (iout,*) "Before scatter npel"
 249       call MPI_Scatter(scount_pmf,1,MPI_INTEGER,npel(it1,it2),1,
 250      &     MPI_INTEGER,Master,WHAM_COMM,ierror)
 251       write (iout,*) "After scatter scount npel",npel(it1,it2)
 252       call flush(iout)
 253       do j=1,4
 254         do i=i1,ii-1
 255           bufin(i)=y(j,i)
 256         enddo
 257         call MPI_Scatterv(bufin(i1),scount_pmf,idispl_pmf,
 258      &    MPI_DOUBLE_PRECISION,
 259      &    bufout(1),npel(it1,it2),MPI_DOUBLE_PRECISION,Master,WHAM_COMM,
 260      &    ierror)
 261         write (iout,*) "After scatter y",j
 262         call flush(iout)
 263         do i=1,npel(it1,it2)
 264           y(j,i+i1-1)=bufout(i)
 265         enddo
 266         do i=i1,ii-1
 267           bufin(i)=w(j,i)
 268         enddo
 269         call MPI_Scatterv(bufin(i1),scount_pmf,idispl_pmf,
 270      &    MPI_DOUBLE_PRECISION,
 271      &    bufout(1),npel(it1,it2),MPI_DOUBLE_PRECISION,Master,WHAM_COMM,
 272      &    ierror)
 273         write (iout,*) "After scatter w",j
 274         call flush(iout)
 275         do i=1,npel(it1,it2)
 276           w(j,i+i1-1)=bufout(i)
 277         enddo
 278       enddo
 279       do j=1,8
 280         do i=i1,ii-1
 281           bufin(i)=zz(j,i)
 282         enddo
 283         call MPI_Scatterv(bufin(i1),scount_pmf,idispl_pmf,
 284      &    MPI_DOUBLE_PRECISION,
 285      &    bufout(1),npel(it1,it2),MPI_DOUBLE_PRECISION,Master,WHAM_COMM,
 286      &    ierror)
 287         do i=1,npel(it1,it2)
 288           zz(j,i+i1-1)=bufout(i)
 289         enddo
 290         write (iout,*) "After scatter z",j
 291         call flush(iout)
 292       enddo
 293       if (lprn) then
 294       write (iout,*) "Distributed eelloc"
 295       do i=i1,i1+npel(it1,it2)-1
 296         write(iout,'(i7,8f7.1,4f10.5)')i,(zz(j,i),j=1,8),(y(j,i),j=1,4)
 297       enddo
 298       endif
 299       ii=i1+npel(it1,it2)
 300       i1=ii
 301 #endif
 302 #endif
 303       write (iout,*)npel(it1,it2)," data points found in file"
 304       write (iout,*)npoint," data points found so far"
 305       call flush(iout)
 306
 307       endif ! eello_pmf
 308
 309       enddo ! it2
 310       enddo ! it1
 311
 312       return
 313    99 write (iout,*) "Error opening PMF fiie(s)"
 314       return1
 315       end
 316 #ifdef MPI
 317 c-------------------------------------------------------------------------------
 318       subroutine work_partition_pmf(n,lprint)
 319       implicit none
 320       include 'DIMENSIONS'
 321       include "DIMENSIONS.ZSCOPT"
 322       include 'mpif.h'
 323       include 'COMMON.MPI'
 324       include 'COMMON.MPISET'
 325       integer ierror,errcode
 326       integer n,ntot,nn,chunk,i,remainder
 327       logical lprint
 328 c      nprocs1=numtasks
 329 C
 330 C Divide conformations between processors; for each proteins the first and
 331 C the last conformation to handle by ith processor is stored in
 332 C indstart(i) and indend(i), respectively.
 333 C
 334 C First try to assign equal number of conformations to each processor.
 335 C
 336       ntot=n
 337       indstart_pmf(0)=1
 338       chunk = N/nprocs1
 339       scount_pmf(0) = chunk
 340       do i=1,nprocs1-1
 341         indstart_pmf(i)=chunk+indstart_pmf(i-1)
 342         scount_pmf(i)=scount_pmf(i-1)
 343       enddo
 344 C
 345 C Determine how many conformations remained yet unassigned.
 346 C
 347       remainder=N-(indstart_pmf(nprocs1-1)+scount_pmf(nprocs1-1)-1)
 348 c      print *,"remainder",remainder
 349 C
 350 C Assign the remainder conformations to consecutive processors, starting
 351 C from the lowest rank; this continues until the list is exhausted.
 352 C
 353       if (remainder .gt. 0) then
 354         do i=1,remainder
 355           scount_pmf(i-1) = scount_pmf(i-1) + 1
 356           indstart_pmf(i) = indstart_pmf(i) + i
 357         enddo
 358         do i=remainder+1,nprocs1-1
 359           indstart_pmf(i) = indstart_pmf(i) + remainder
 360         enddo
 361       endif
 362
 363       indstart_pmf(nprocs1)=N+1
 364       scount_pmf(nprocs1)=0
 365
 366       do i=0,NProcs1
 367         indend_pmf(i)=indstart_pmf(i)+scount_pmf(i)-1
 368         idispl_pmf(i)=indstart_pmf(i)-1
 369       enddo
 370
 371       N=0
 372       do i=0,Nprocs1-1
 373         N=N+indend_pmf(i)-indstart_pmf(i)+1
 374       enddo
 375
 376       if (N.ne.ntot) then
 377         write (2,*) "!!! Checksum error on processor",me
 378         call flush(2)
 379         call MPI_Abort( MPI_COMM_WORLD, Ierror, Errcode )
 380       endif
 381
 382       if (lprint) then
 383         write (2,*) "Partition of work between processors"
 384         do i=0,nprocs1-1
 385           write (2,'(a,i5,a,i7,a,i7,a,i7)')
 386      &      "Processor",i," indstart",indstart_pmf(i),
 387      &      " indend",indend_pmf(i)," count",scount_pmf(i)," idispl",
 388      &      idispl_pmf(i)
 389         enddo
 390       endif
 391       return
 392       end
 393 #endif