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