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