source/maxlik/src-Fmatch/proc.F

   1       subroutine proc_data(nvarr,x,*)
   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),tag
   8       include "COMMON.MPI"
   9 #endif
  10       include "COMMON.CONTROL"
  11       include "COMMON.ENERGIES"
  12       include "COMMON.VMCPAR"
  13       include "COMMON.OPTIM"
  14       include "COMMON.WEIGHTS"
  15       include "COMMON.PROTNAME"
  16       include "COMMON.IOUNITS"
  17       include "COMMON.FFIELD"
  18       include "COMMON.PROTFILES"
  19       include "COMMON.COMPAR"
  20       include "COMMON.CHAIN"
  21       include "COMMON.CLASSES"
  22       include "COMMON.THERMAL"
  23       include "COMMON.WEIGHTDER"
  24       include "COMMON.EREF"
  25       include "COMMON.SBRIDGE"
  26       include "COMMON.INTERACT"
  27       include "COMMON.VAR"
  28       include "COMMON.TIME1"
  29       include "COMMON.FMATCH"
  30       double precision fac0,t0
  31       double precision fT(5),ftprim(5),ftbis(5),quot,quotl,quotl1,
  32      &  denom,kfac,kfacl,logfac,eplus,eminus,tanhT,entfacmax,
  33      &  sument,sumentsq,aux
  34       real csingle(3,maxres2)
  35       double precision creff(3,maxres2),cc(3,maxres2)
  36       double precision thetareff(maxres2),phireff(maxres2),
  37      & xxreff(maxres2),yyreff(maxres2),zzreff(maxres2)
  38       double precision przes(3),obrot(3,3),etoti
  39       logical non_conv
  40       integer wykl
  41       integer iprot,i,ind,ii,j,k,kk,l,ncl,is,ie,nsum,nng,nvarr,ng,
  42      &  inn,in,iref,ib,ibatch,num,itmp,imin
  43       integer l1,l2,ncon_all
  44       double precision rcyfr
  45       integer jj,istart_conf,iend_conf,ipass_conf,iii
  46       integer ilen,iroof
  47       external ilen,iroof
  48       double precision ej,emaxj,sumP,sumP_t
  49       double precision fac
  50       double precision x(max_paropt)
  51       double precision wconf(maxstr_proc)
  52       double precision dx,dy,dz,dxref,dyref,dzref,
  53      &  dtheta,dphi,rmsthet,rmsphi,rmsside
  54       double precision Fmean,SStot
  55       double precision pinorm
  56       external pinorm
  57 #ifdef MPI
  58       double precision e_total_T_(maxstr_proc)
  59 #endif
  60       logical LPRN
  61       integer maxl1,maxl2
  62       double precision rms,rmsave,rmsmin,rmsnat,qwolynes,qwolyness
  63       cutoffeval=.false.
  64       T0=3.0d2
  65       fac0=1.0d0/(T0*Rgas)
  66       kfac=2.4d0
  67       write (iout,*) "fac0",fac0
  68       write(iout,*) "torsional and valence angle mode",tor_mode
  69 c
  70 c Determine the number of variables and the initial vector of optimizable
  71 c parameters
  72 c
  73       lprn=.false.
  74       if (restart_flag) then
  75         call w2x(nvarr,x_orig,*11)
  76         do i=1,n_ene
  77           ww_oorig(i)=ww(i)
  78         enddo
  79         open (88,file=prefix(:ilen(prefix))//'.restin',status='unknown')
  80         read(88,*,end=99,err=99) (x(i),i=1,nvarr)
  81         close(88)
  82         write (iout,*)
  83         write (iout,*)
  84      &    "Variables replaced with values from the restart file."
  85         do i=1,nvarr
  86           write (iout,'(i5,f15.5)') i,x(i)
  87         enddo
  88         write (iout,*)
  89         call x2w(nvarr,x)
  90         goto 88
  91   99    write (iout,*)
  92      &    "Error - restart file nonexistent, incompatible or damaged."
  93 #ifdef MPI
  94         call MPI_Finalize(Ierror)
  95 #endif
  96         stop
  97      &    "Error - restart file nonexistent, incompatible or damaged."
  98       else
  99         call w2x(nvarr,x,*11)
 100         write (iout,*) "PROC: after W2X nvarr",nvarr
 101         do i=1,n_ene
 102           ww_oorig(i)=ww(i)
 103         enddo
 104         do i=1,nvarr
 105           x_orig(i)=x(i)
 106         enddo
 107       endif
 108   88  continue
 109 c
 110 c Setup weights for UNRES
 111 c
 112       wsc=ww(1)
 113       wscp=ww(2)
 114       welec=ww(3)
 115       wcorr=ww(4)
 116       wcorr5=ww(5)
 117       wcorr6=ww(6)
 118       wel_loc=ww(7)
 119       wturn3=ww(8)
 120       wturn4=ww(9)
 121       wturn6=ww(10)
 122       wang=ww(11)
 123       wscloc=ww(12)
 124       wtor=ww(13)
 125       wtor_d=ww(14)
 126       wstrain=ww(15)
 127       wvdwpp=ww(16)
 128       wbond=ww(17)
 129       wsccor=ww(19)
 130 #ifdef SCP14
 131       scal14=ww(17)/ww(2)
 132 #endif
 133       do i=1,n_ene
 134         weights(i)=ww(i)
 135       enddo
 136
 137       DO IPROT=1,NPROT
 138
 139       if (.not.maxlik(iprot)) cycle
 140
 141       call restore_molinfo(iprot)
 142
 143
 144       DO IBATCH=1,natlike(iprot)+2
 145
 146 c Calculate temperature-dependent weight scaling factors
 147       do ib=1,nbeta(ibatch,iprot)
 148         fac=betaT(ib,ibatch,iprot)
 149         quot=fac0/fac
 150         if (rescale_mode.eq.0) then
 151           do l=1,5
 152             fT(l)=1.0d0
 153             fTprim(l)=0.0d0
 154             fTbis(l)=0.0d0
 155           enddo
 156         else if (rescale_mode.eq.1) then
 157           quotl=1.0d0
 158           kfacl=1.0d0
 159           do l=1,5
 160             quotl1=quotl
 161             quotl=quotl*quot
 162             kfacl=kfacl*kfac
 163             denom=kfacl-1.0d0+quotl
 164             fT(l)=kfacl/denom
 165             ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 166             ftbis(l)=l*kfacl*quotl1*
 167      &       (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 168             write (iout,*) "ib",ib," beta",betaT(ib,ibatch,iprot),
 169      &      " quot",quot," l",l,fT(l),fTprim(l),fTbis(l)
 170           enddo
 171         else if (rescale_mode.eq.2) then
 172           quotl=1.0d0
 173           do l=1,5
 174             quotl1=quotl
 175             quotl=quotl*quot
 176             eplus=dexp(quotl)
 177             eminus=dexp(-quotl)
 178             logfac=1.0d0/dlog(eplus+eminus)
 179             tanhT=(eplus-eminus)/(eplus+eminus)
 180             fT(l)=1.12692801104297249644d0*logfac
 181             ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
 182             ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
 183      &      2*l*quotl1/T0*logfac*
 184      &      (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
 185      &      +ftprim(l)*tanhT)
 186           enddo
 187         else
 188           write (iout,*) "Unknown RESCALE_MODE",rescale_mode
 189           call flush(iout)
 190           stop
 191         endif
 192         do k=1,n_ene
 193           escal(k,ib,ibatch,iprot)=1.0d0
 194           escalprim(k,ib,ibatch,iprot)=0.0d0
 195           escalbis(k,ib,ibatch,iprot)=0.0d0
 196         enddo
 197         if (shield_mode.gt.0) then
 198           escal(1,ib,ibatch,iprot)=ft(1)
 199           escal(2,ib,ibatch,iprot)=ft(1)
 200           escal(16,ib,ibatch,iprot)=ft(1)
 201           escalprim(1,ib,ibatch,iprot)=ft(1)
 202           escalprim(2,ib,ibatch,iprot)=ft(1)
 203           escalprim(16,ib,ibatch,iprot)=ft(1)
 204           escalbis(1,ib,ibatch,iprot)=ft(1)
 205           escalbis(2,ib,ibatch,iprot)=ft(1)
 206           escalbis(16,ib,ibatch,iprot)=ft(1)
 207         endif
 208         escal(3,ib,ibatch,iprot)=ft(1)
 209         escal(4,ib,ibatch,iprot)=ft(3)
 210         escal(5,ib,ibatch,iprot)=ft(4)
 211         escal(6,ib,ibatch,iprot)=ft(5)
 212         escal(7,ib,ibatch,iprot)=ft(2)
 213         escal(8,ib,ibatch,iprot)=ft(2)
 214         escal(9,ib,ibatch,iprot)=ft(3)
 215         escal(10,ib,ibatch,iprot)=ft(5)
 216         escal(13,ib,ibatch,iprot)=ft(1)
 217         escal(14,ib,ibatch,iprot)=ft(2)
 218         escal(19,ib,ibatch,iprot)=ft(1)
 219         escalprim(3,ib,ibatch,iprot)=ftprim(1)
 220         escalprim(4,ib,ibatch,iprot)=ftprim(3)
 221         escalprim(5,ib,ibatch,iprot)=ftprim(4)
 222         escalprim(6,ib,ibatch,iprot)=ftprim(5)
 223         escalprim(7,ib,ibatch,iprot)=ftprim(2)
 224         escalprim(8,ib,ibatch,iprot)=ftprim(2)
 225         escalprim(9,ib,ibatch,iprot)=ftprim(3)
 226         escalprim(10,ib,ibatch,iprot)=ftprim(5)
 227         escalprim(13,ib,ibatch,iprot)=ftprim(1)
 228         escalprim(14,ib,ibatch,iprot)=ftprim(2)
 229         escalprim(19,ib,ibatch,iprot)=ftprim(1)
 230         escalbis(3,ib,ibatch,iprot)=ftbis(1)
 231         escalbis(4,ib,ibatch,iprot)=ftbis(3)
 232         escalbis(5,ib,ibatch,iprot)=ftbis(4)
 233         escalbis(6,ib,ibatch,iprot)=ftbis(5)
 234         escalbis(7,ib,ibatch,iprot)=ftbis(2)
 235         escalbis(8,ib,ibatch,iprot)=ftbis(2)
 236         escalbis(9,ib,ibatch,iprot)=ftbis(3)
 237         escalbis(10,ib,ibatch,iprot)=ftbis(5)
 238         escalbis(13,ib,ibatch,iprot)=ftbis(1)
 239         escalbis(14,ib,ibatch,iprot)=ftbis(2)
 240         escalbis(19,ib,ibatch,iprot)=ftbis(1)
 241       enddo
 242
 243       betmin(ibatch,iprot)=betaT(1,ibatch,iprot)
 244       betmax(ibatch,iprot)=betaT(1,ibatch,iprot)
 245       do ib=2,nbeta(ibatch,iprot)
 246         if (betaT(ib,ibatch,iprot).lt.betmin(ibatch,iprot))
 247      &     betmin(ibatch,iprot)=betaT(ib,ibatch,iprot)
 248         if (betaT(ib,ibatch,iprot).gt.betmax(ibatch,iprot))
 249      &     betmax(ibatch,iprot)=betaT(ib,ibatch,iprot)
 250         enddo
 251        write (iout,'(2(a,i5,2x),2(a,f7.2,2x))')
 252      &    "Protein",iprot," batch",ibatch,
 253      &    " betmin",betmin(ibatch,iprot)," betmax",betmax(ibatch,iprot)
 254
 255       ENDDO ! IBATCH
 256
 257       ENDDO ! IPROT
 258 c
 259 c Compute energy and entropy histograms to filter out conformations
 260 c with potntially too low or too high weights.
 261 c
 262       call determine_histograms
 263 c
 264 c Make the working list of conformations
 265 c
 266       call make_list(.true.,.true.,nvarr,x)
 267
 268       call make_list_MD(.true.)
 269
 270 c 3/8/2013 Calculate RMSDs/Qs and the Ptab array
 271       call PTab_calc
 272 c AL 5/13/2019 Compute variances of MD forces
 273       do iprot=1,nprot
 274
 275         if (.not.fmatch(iprot)) cycle
 276
 277         open (icbase,file=bprotfiles_MD(iprot),status="old",
 278      &    form="unformatted",access="direct",recl=lenrec_MD(iprot))
 279
 280         Fmean=0.0d0
 281         do ii=1,ntot_work_MD(iprot),maxstr_proc
 282           istart_conf=ii
 283           iend_conf=min0(ii+maxstr_proc-1,ntot_MD(iprot))
 284           call daread_MDcoords(iprot,istart_conf,iend_conf)
 285           do j=istart_conf,iend_conf
 286 c            jj=list_conf_MD(j,iprot)
 287 c          write (iout,*) "iprot",iprot," j",j," jj",jj
 288             do i=1,nsite(iprot)
 289               do k=1,3
 290                 Fmean=Fmean+CGforce_MD(k,i,j,iprot)
 291               enddo
 292             enddo
 293           enddo
 294         enddo
 295
 296         Fmean=Fmean/(3*nsite(iprot)*ntot_work_MD(iprot))
 297         write(iout,*) "PROC_DATA Protein",iprot," Fmean",Fmean
 298
 299         SStot=0.0d0
 300
 301         do j=1,ntot_work_MD(iprot)
 302 c          jj=list_conf_MD(j,iprot)
 303           do i=1,nsite(iprot)
 304             do k=1,3
 305               SStot=SStot+(CGforce_MD(k,i,j,iprot)-Fmean)**2
 306             enddo
 307           enddo
 308         enddo
 309
 310         varforce(iprot)=SStot/(3*nsite(iprot)*ntot_work_MD(iprot))
 311         meanforce(iprot)=Fmean
 312         write (iout,*) "PROC_DATA Protein",iprot," varforce",
 313      &   varforce(iprot)," mean force",meanforce(iprot)
 314
 315         close(icbase)
 316
 317       enddo ! iprot
 318
 319       return
 320    11 return1
 321       end
 322 c-------------------------------------------------------------------------------
 323       subroutine determine_histograms
 324       implicit none
 325       include "DIMENSIONS"
 326       include "DIMENSIONS.ZSCOPT"
 327       integer maxind
 328       parameter (maxind=1000)
 329 #ifdef MPI
 330       include "mpif.h"
 331       integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
 332       include "COMMON.MPI"
 333 #endif
 334       include "COMMON.COMPAR"
 335       include "COMMON.ENERGIES"
 336       include "COMMON.VMCPAR"
 337       include "COMMON.WEIGHTS"
 338       include "COMMON.PROTNAME"
 339       include "COMMON.IOUNITS"
 340       include "COMMON.FFIELD"
 341       include "COMMON.PROTFILES"
 342       include "COMMON.CHAIN"
 343       include "COMMON.CONTROL"
 344       include "COMMON.CLASSES"
 345       integer iprot,i,ii,iii,j,k,l,ibatch,inde,indt,idelte,ideltt
 346       double precision eminh,emaxh,tminh,tmaxh
 347       integer he(0:maxind),
 348      &  ht(0:maxind),hemax,htmax
 349 c
 350 c Construct energy and entropy histograms
 351 c
 352       write (iout,*) "Calling DETERMINE HISTOGRAMS"
 353       call flush(iout)
 354       do iprot=1,nprot
 355           if (.not.maxlik(iprot)) cycle
 356           eminh=1.0d20
 357           emaxh=-1.0d20
 358           tminh=1.0d20
 359           tmaxh=-1.0d20
 360           do i=0,maxind
 361             he(i)=0
 362             ht(i)=0
 363           enddo
 364           do i=1,ntot(iprot)
 365             iii=i
 366             if (eini(iii,iprot).lt.eminh)
 367      &        eminh=eini(iii,iprot)
 368             if (eini(iii,iprot).gt.emaxh)
 369      &        emaxh=eini(iii,iprot)
 370              if (entfac(iii,iprot).lt.tminh)
 371      &        tminh=entfac(iii,iprot)
 372             if (entfac(iii,iprot).gt.tmaxh)
 373      &        tmaxh=entfac(iii,iprot)
 374           enddo
 375 c          write (2,*)  "DETERMINE HISTOGRAMS 1"
 376           call flush(iout)
 377           do i=1,ntot(iprot)
 378             iii=i
 379             inde=eini(iii,iprot)-eminh
 380             indt=entfac(iii,iprot)-tminh
 381             if (inde.le.maxind) he(inde)=he(inde)+1
 382             if (indt.le.maxind) ht(indt)=ht(indt)+1
 383           enddo
 384 c          write (2,*)  "DETERMINE HISTOGRAMS 2"
 385             idelte=emaxh-eminh
 386             ideltt=tmaxh-tminh
 387             hemax=0
 388 c            write (iout,*) "idelte",idelte," ideltt",ideltt
 389             call flush(iout)
 390             do i=0,idelte
 391 c              write (iout,*) "i",i," he",he(i)," hemax",hemax
 392               call flush(iout)
 393               if (he(i).gt.hemax) hemax=he(i)
 394             enddo
 395             htmax=0
 396             do i=0,ideltt
 397 c              write (iout,*) "i",i," ht",ht(i)," htmax",htmax
 398               call flush(iout)
 399               if (ht(i).gt.htmax) htmax=ht(i)
 400             enddo
 401 c          write (2,*)  "DETERMINE HISTOGRAMS 3"
 402           call flush(iout)
 403             hemax=hemax/hefac(iprot)
 404             if (hemax.lt.hemax_low(iprot))
 405      &        hemax=hemax_low(iprot)
 406             htmax=htmax/htfac(iprot)
 407             if (htmax.lt.htmax_low(iprot))
 408      &        htmax=htmax_low(iprot)
 409             i=0
 410             do while (i.lt.idelte .and. he(i).lt.hemax)
 411               i=i+1
 412             enddo
 413 c          write (2,*)  "DETERMINE HISTOGRAMS 4"
 414           call flush(iout)
 415             e_lowb(iprot)=eminh+i
 416             i=0
 417             do while (i.lt.ideltt .and. ht(i).lt.htmax)
 418               i=i+1
 419             enddo
 420             t_lowb(iprot)=tminh+i
 421 #ifdef DEBUG
 422             write (iout,*) "protein",iprot
 423             write (iout,*) "energy histogram"
 424             do i=0,idelte
 425               write (iout,'(f15.5,i10)') eminh+i,he(i)
 426             enddo
 427             write (iout,*) "entropy histogram"
 428             do i=0,ideltt
 429               write (iout,'(f15.5,i10)') tminh+i,ht(i)
 430             enddo
 431             write (iout,*) "emin",eminh," tmin",tminh,
 432      &       " hemax",hemax," htmax",htmax,
 433      &       " e_lowb",e_lowb(iprot),
 434      &       " t_lowb",t_lowb(iprot)
 435             call flush(iout)
 436 #endif
 437       enddo
 438       return
 439       end
 440 c------------------------------------------------------------------------------
 441       subroutine imysort(n, x, ipermut)
 442       implicit none
 443       integer n
 444       integer x(n),xtemp
 445       integer ipermut(n)
 446       integer i,j,imax,ipm
 447       do i=1,n
 448         xtemp=x(i)
 449         imax=i
 450         do j=i+1,n
 451           if (x(j).lt.xtemp) then
 452             imax=j
 453             xtemp=x(j)
 454           endif
 455         enddo
 456         x(imax)=x(i)
 457         x(i)=xtemp
 458         ipm=ipermut(imax)
 459         ipermut(imax)=ipermut(i)
 460         ipermut(i)=ipm
 461       enddo
 462       return
 463       end
 464 c--------------------------------------------------------------------------
 465       subroutine write_conf_count
 466       implicit none
 467       include "DIMENSIONS"
 468       include "DIMENSIONS.ZSCOPT"
 469       include "COMMON.CLASSES"
 470       include "COMMON.COMPAR"
 471       include "COMMON.VMCPAR"
 472       include "COMMON.PROTNAME"
 473       include "COMMON.IOUNITS"
 474       include "COMMON.PROTFILES"
 475       integer iprot,ibatch,i,ii,j,kk
 476       integer ilen
 477       external ilen
 478       logical LPRN
 479       write (iout,*)
 480       write (iout,*)"Numbers of conformations after applying the cutoff"
 481       write (iout,*)
 482       do iprot=1,nprot
 483         write (iout,'(a,2x,a,i10)') "Protein",
 484      &    protname(iprot)(:ilen(protname(iprot))),
 485      &    ntot_work(iprot)
 486       enddo
 487       return
 488       end
 489 c-------------------------------------------------------------------------
 490       double precision function fdum()
 491       fdum=0.0D0
 492       return
 493       end