source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/func1_sc.F

   1       subroutine func1(n,What,x)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       integer IERROR,ErrCode
   8       include "COMMON.MPI"
   9 #endif
  10       integer n,nf,What
  11       double precision x(n)
  12       include "COMMON.WEIGHTS"
  13       include "COMMON.WEIGHTDER"
  14       include "COMMON.CLASSES"
  15       include "COMMON.ENERGIES"
  16       include "COMMON.IOUNITS"
  17       include "COMMON.VMCPAR"
  18       include "COMMON.OPTIM"
  19       include "COMMON.NAMES"
  20       include "COMMON.INTERACT"
  21       include "COMMON.TIME1"
  22       include "COMMON.CHAIN"
  23       include "COMMON.PROTFILES"
  24       include "COMMON.COMPAR"
  25       include "COMMON.PROTNAME"
  26       include "COMMON.THERMAL"
  27       include "COMMON.TORSION"
  28 C Define local variables
  29       integer i,ii,iii,ind,jj,j,k,l,ik,jk,iprot,ib,l1,l2,ll1,ll2,
  30      &  it
  31       integer ibatch
  32       integer ipass_conf,istart_conf,iend_conf
  33       double precision energia(0:max_ene)
  34       double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk,
  35      &  eave_pft(max_ene),emix_pft(max_ene),
  36      &  elowest_t(2,maxT),esquare_ft,
  37      &  elowest_aux(2,maxT),
  38      &  efree_t,emixsq_pft(max_ene),
  39      &  eneps_mixt(nntyp),eneps_meant(nntyp),
  40      &  enepsave_ft(nntyp),eneps_mix_ft(nntyp),
  41      &  eneps_mixsq_ft(nntyp),emean_ft,fac
  42 #ifdef MPI
  43       double precision e_total_(maxstr_proc),eini_(maxstr_proc),
  44      &  entfac_(maxstr_proc),rmstb_(maxstr_proc)
  45 #endif
  46       double precision aux
  47       double precision tcpu_ini,tcpu
  48       logical lprn
  49       integer icant
  50       external icant
  51       character*80 vormat
  52
  53       tcpu_ini = tcpu()
  54       n_func=n_func+1
  55       lprn=.true.
  56 c#define DEBUG
  57 #ifdef DEBUG
  58       write (iout,*) "FUNC1: Init_Ene ",init_ene
  59       call flush(iout)
  60       write (iout,*) "ELOWEST at the beginning of FUC1"
  61       do iprot=1,nprot
  62         do ibatch=1,natlike(iprot)+2
  63           do ib=1,nbeta(ibatch,iprot)
  64             write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
  65      &         " elowest",elowest(ib,ibatch,iprot)
  66           enddo
  67         enddo
  68       enddo
  69 #endif
  70 #ifdef MPI
  71       call MPI_Bcast( What, 1, MPI_INTEGER, Master, Comm1, IERROR)
  72       call MPI_Bcast( cutoffeval, 1, MPI_LOGICAL, Master, Comm1, IERROR)
  73       call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION, Master, Comm1, IERROR)
  74       if (What.le.0) return
  75 #ifdef DEBUG
  76         write (iout,*) "Processor",me,me1," What",What
  77         write (iout,*) "Processor",me,me1," x",(x(i),i=1,n)
  78         call flush(iout)
  79 #endif
  80 #else
  81 #ifdef DEBUG
  82         write (iout,*) "x",(x(i),i=1,n)
  83 #endif
  84 #endif
  85 c#undef DEBUG
  86 C Convert variables into weights and other UNRES energy parameters.
  87       call x2w(n,x)
  88 C Calculate the total energies.
  89
  90 c      write (iout,*) "Processor",me," nprot",nprot
  91 c Revise the list of conformations
  92       if (What.eq.5) then
  93         write (iout,*) "MAKE_LIST called"
  94         call flush(iout)
  95         call MPI_Bcast( enecut(1), nprot, MPI_DOUBLE_PRECISION,
  96      &    Master, Comm1, IERROR)
  97         call make_list(.true.,.false.,n,x)
  98         return
  99       endif
 100       DO iprot=1,nprot
 101 #ifdef DEBUG
 102         if (What.eq.2) then
 103           open(icsa_rbank,file="corr."//protname(iprot),
 104      &    access="direct",form="formatted",recl=25+10*natlike(iprot))
 105           vormat='(2i6,2i3,e15.5,'
 106           if (natlike(iprot).lt.10) then
 107            write(vormat(16:),'(i1,a)') maxdimnat*natlike(iprot),'f7.4/)'
 108           else
 109            write(vormat(16:),'(i2,a)') maxdimnat*natlike(iprot),'f7.4/)'
 110           endif
 111           write (iout,*) "format ",vormat
 112           call flush(iout)
 113         endif
 114 #endif
 115 #ifdef DEBUG
 116       write (iout,*) "Processor",me," iprot",iprot,
 117      & " indstart",indstart(me1,iprot)," indend",indend(me1,iprot),
 118      & " init_ene",init_ene," mod_fourier",mod_fourier(nloctyp),
 119      & " mod_fouriertor",mod_fouriertor(nloctyp),
 120      & " mod_elec",mod_elec," mod_scp",mod_scp," mod_tor",mod_tor,
 121      & " mod_sccor",mod_sccor," mod_ang",mod_ang
 122 #endif
 123       call restore_molinfo(iprot)
 124       ii=0
 125
 126       IF (NCHUNK_CONF(IPROT).EQ.1) THEN
 127
 128 #ifdef MPI
 129       do i=indstart(me1,iprot),indend(me1,iprot)
 130 #else
 131       do i=1,ntot_work(iprot)
 132 #endif
 133 c        write (iout,*) "i",i," ii",ii," calling ene_eval"
 134         call ene_eval(iprot,i,ii)
 135 c        write (iout,*) "After ene_eval: i",i," ii",ii
 136 #ifdef DEBUG
 137         if (What.eq.2) then
 138           write (icsa_rbank,vormat,rec=i)
 139      &     i,list_conf(i,iprot),
 140      &    e_total(i,iprot),((nu(l,k,ii,iprot),l=1,natdim(k,iprot)),
 141      &    k=1,natlike(iprot))
 142         endif
 143 #endif
 144       enddo
 145
 146       ELSE
 147
 148       if (.not.init_ene .and. mod_fourier(nloctyp).eq.0
 149      &  .and. mod_fouriertor(nloctyp).eq.0
 150      &  .and. .not. mod_elec .and. .not. mod_scp
 151      &  .and. .not. mod_tor .and. .not. mod_sccor) then
 152         open (ientin,file=benefiles(iprot),status="old",
 153      &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
 154         if (natlike(iprot).gt.0)
 155      &    open (icbase,file=bprotfiles(iprot),status="old",
 156      &    form="unformatted",access="direct",recl=lenrec(iprot))
 157       else
 158         open (icbase,file=bprotfiles(iprot),status="old",
 159      &    form="unformatted",access="direct",recl=lenrec(iprot))
 160 c        write (iout,*) "lenrec_ene",lenrec_ene(iprot)
 161         call flush(iout)
 162 c Change AL 12/30/2017
 163         if (.not. mod_other_params)
 164      &  open (ientout,file=benefiles(iprot),status="unknown",
 165      &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
 166       endif
 167       ipass_conf=0
 168 #ifdef MPI
 169       do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
 170         istart_conf=i
 171         iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
 172 #else
 173       do i=1,ntot_work(iprot),maxstr_proc
 174         istart_conf=i
 175         iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
 176 #endif
 177         ii=0
 178 c
 179 c Read the chunk of conformations off a DA scratchfile.
 180 c
 181         if (.not. init_ene .and. mod_fourier(nloctyp).eq.0
 182      &  .and. mod_fouriertor(nloctyp).eq.0
 183      &  .and. .not. mod_elec .and. .not. mod_scp
 184      &  .and. .not. mod_tor .and. .not. mod_sccor .and. .not. mod_ang)
 185      &  then
 186 #ifdef DEBUG
 187           write (iout,*) "func1: calling daread_ene",istart_conf,
 188      &      iend_conf
 189 #endif
 190           call daread_ene(iprot,istart_conf,iend_conf)
 191           if (natlike(iprot).gt.0) then
 192 #ifdef DEBUG
 193             write (iout,*) "Calling daread_coords",istart_conf,
 194      &       iend_conf
 195 #endif
 196             call daread_ccoords(iprot,istart_conf,iend_conf)
 197           endif
 198           do iii=istart_conf,iend_conf
 199             call ene_eval(iprot,iii,ii)
 200           enddo
 201 #ifdef DEBUG
 202           write (iout,*) "func1: calling dawrite_ene",istart_conf,
 203      &      iend_conf
 204 #endif
 205           call dawrite_ene(iprot,istart_conf,iend_conf,ientin)
 206         else
 207           call daread_ccoords(iprot,istart_conf,iend_conf)
 208           do iii=istart_conf,iend_conf
 209             call ene_eval(iprot,iii,ii)
 210           enddo
 211 c AL 12/30/17 - writing energies is not needed when whole energy is evaluated
 212 c          call dawrite_ene(iprot,istart_conf,iend_conf,ientout)
 213         endif
 214 #ifdef DEBUG
 215         if (What.eq.2) then
 216           do iii=istart_conf,iend_conf
 217             write (icsa_rbank,vormat,rec=iii)
 218      &      iii,list_conf(iii,iprot),
 219      &      e_total(iii,iprot),
 220      &      ((nu(l,k,iii-istart_conf+1,iprot),l=1,natdim(k,iprot)),
 221      &      k=1,natlike(iprot))
 222           enddo
 223         endif
 224 #endif
 225       enddo
 226
 227       if (.not.init_ene .and. mod_fourier(nloctyp).eq.0
 228      &  .and. mod_fouriertor(nloctyp).eq.0
 229      &  .and. .not. mod_elec .and. .not. mod_scp
 230      &  .and. .not. mod_tor .and. .not. mod_sccor .and. .not. mod_ang)
 231      &   then
 232         close (ientin)
 233         if (natlike(iprot).gt.0) close(icbase)
 234       else
 235         close (icbase)
 236       endif
 237       if (init_ene .or. mod_fourier(nloctyp).gt.0 .or. mod_elec
 238      &   .or. mod_fouriertor(nloctyp).gt.0
 239      &   .or. mod_scp .or. mod_tor .or. mod_sccor .or. mod_ang)
 240      &  close(ientout)
 241
 242       ENDIF
 243 #ifdef MPI
 244 c      print *,"Processor",me,me1," finished computing energies"
 245 c      call flush(iout)
 246       do i=1,scount(me1,iprot)
 247         e_total_(i)=e_total(indstart(me1,iprot)+i-1,iprot)
 248         eini_(i)=eini(indstart(me1,iprot)+i-1,iprot)
 249         entfac_(i)=entfac(indstart(me1,iprot)+i-1,iprot)
 250         rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot)
 251       enddo
 252       if (What.eq.2 .or. What.eq.3) then
 253 c        call MPI_Gatherv(e_total(indstart(me1,iprot),iprot),
 254 c     &     scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
 255 c     &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 256 c     &     Comm1, IERROR)
 257         call MPI_Gatherv(e_total_(1),
 258      &     scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
 259      &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 260      &     Comm1, IERROR)
 261 c        call MPI_Gatherv(entfac(indstart(me1,iprot),iprot),
 262 c     &     scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot),
 263 c     &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 264 c     &     Comm1, IERROR)
 265         call MPI_Gatherv(entfac_(1),
 266      &     scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot),
 267      &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 268      &     Comm1, IERROR)
 269 c        call MPI_Gatherv(eini(indstart(me1,iprot),iprot),
 270 c     &     scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot),
 271 c     &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 272 cc     &     Comm1, IERROR)
 273         call MPI_Gatherv(eini_(1),
 274      &     scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot),
 275      &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 276      &     Comm1, IERROR)
 277 c        call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
 278 c     &     scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
 279 c     &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 280 c     &     Comm1, IERROR)
 281         call MPI_Gatherv(rmstb_(1),
 282      &     scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
 283      &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 284      &     Comm1, IERROR)
 285       endif
 286       if (What.eq.3) then
 287 c         call MPI_Allgatherv(e_total(indstart(me1,iprot),iprot),
 288 c     &     scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
 289 c     &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
 290 c     &     Comm1, IERROR)
 291          call MPI_Allgatherv(e_total_(1),
 292      &     scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
 293      &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
 294      &     Comm1, IERROR)
 295       endif
 296 #endif
 297
 298 c Mean and lowest energies of structural classes
 299
 300       IF (NCHUNK_CONF(IPROT).EQ.1) THEN
 301
 302 #ifdef MPI
 303       do i=1,ntot_work(iprot)
 304         i2ii(i,iprot)=0
 305       enddo
 306       ii=0
 307       do i=indstart(me1,iprot),indend(me1,iprot)
 308         ii=ii+1
 309         i2ii(i,iprot)=ii
 310       enddo
 311       istart_conf=indstart(me1,iprot)
 312       iend_conf=indend(me1,iprot)
 313 #else
 314       do i=1,ntot_work(iprot)
 315         i2ii(i,iprot)=i
 316       endif
 317       istart_conf=1
 318       iend_conf=ntot_work(iprot)
 319 #endif
 320 #ifdef DEBUG
 321       do i=1,ntot_work(iprot)
 322         write (iout,*) "i",i," i2ii",i2ii(i,iprot)
 323       enddo
 324       call flush(iout)
 325 #endif
 326 C Change AL 12/30/2017
 327       if (.not.mod_other_params)
 328      &open (ientin,file=benefiles(iprot),status="old",
 329      &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
 330 c      call daread_ene(iprot,istart_conf,iend_conf)
 331       call emin_search(iprot)
 332
 333       ELSE
 334
 335       if (.not.mod_other_params)
 336      &open (ientin,file=benefiles(iprot),status="old",
 337      &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
 338       ipass_conf=0
 339 #ifdef MPI
 340       do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
 341         iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
 342 #else
 343       do istart_conf=1,ntot_work(iprot),maxstr_proc
 344         iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
 345 #endif
 346 c
 347 c Read the chunk of energies and derivatives off a DA scratchfile.
 348 c
 349         ipass_conf=ipass_conf+1
 350         do i=1,ntot_work(iprot)
 351           i2ii(i,iprot)=0
 352         enddo
 353         ii=0
 354         do i=istart_conf,iend_conf
 355           ii=ii+1
 356           i2ii(i,iprot)=ii
 357         enddo
 358 #ifdef DEBUG
 359         write (iout,*) "ipass_conf",ipass_conf,
 360      &    " istart_conf",istart_conf," iend_conf",iend_conf
 361         do i=1,ntot_work(iprot)
 362           write (iout,*) "i",i," i2ii",i2ii(i,iprot)
 363         enddo
 364         call flush(iout)
 365 #endif
 366         call daread_ene(iprot,istart_conf,iend_conf)
 367         call emin_search(iprot)
 368       enddo
 369
 370       close(ientin)
 371       ENDIF
 372
 373       ENDDO ! iprot
 374
 375 #ifdef MPI
 376 c Complete the calculation of the lowest energies over all classes and
 377 c distribute the values to all procs
 378       do iprot=1,nprot
 379 #ifdef DEBUG
 380         write (iout,*) "Processor",me,me1," calling MPI_Reduce: 2"
 381         do k=1,natlike(iprot)+2
 382         write (iout,*) "iprot",iprot," k",k
 383         do ib=1,nbeta(k,iprot)
 384         write (iout,*) "ib",ib
 385         write (iout,'(7hELOWEST,10e15.3)')
 386      &    elowest(ib,k,iprot)
 387         enddo
 388         enddo
 389         call flush(iout)
 390 #endif
 391         do ibatch=1,natlike(iprot)+2
 392           do ib=1,nbeta(ibatch,iprot)
 393             elowest_aux(1,ib)=elowest(ib,ibatch,iprot)
 394             elowest_aux(2,ib)=ind_lowest(ib,ibatch,iprot)
 395           enddo
 396           call MPI_Allreduce(elowest_aux(1,1),elowest_t(1,1),
 397      &      nbeta(ibatch,iprot),
 398      &      MPI_2DOUBLE_PRECISION, MPI_MINLOC, Comm1, IERROR)
 399 #ifdef DEBUG
 400           do ib=1,nbeta(ibatch,iprot)
 401             write (iout,*) "beta=",betaT(ib,ibatch,iprot)
 402             write (iout,'(9helowest_t,10f15.3)')
 403      &        elowest_t(1,ib),elowest_t(2,ib)
 404           enddo
 405           write (iout,*) "Processor",me,me1," finished MPI_Reduce: 2"
 406 #endif
 407           do ib=1,nbeta(ibatch,iprot)
 408 c            write (iout,*) "IB",ib," ibatch",ibatch," iprot",iprot
 409             call flush(iout)
 410             elowest(ib,ibatch,iprot)=elowest_t(1,ib)
 411 c            write (iout,*) "elowest",ib,ibatch,iprot,
 412 c     &              elowest(ib,ibatch,iprot)
 413             call flush(iout)
 414             ind_lowest(ib,ibatch,iprot)=elowest_t(2,ib)
 415           enddo
 416         enddo ! ibatc
 417
 418       enddo ! iprot
 419 #endif
 420       do iprot=1,nprot
 421         call averages(iprot)
 422 c        write (iout,*) "After averages, calling Thermal"
 423 c        call flush(iout)
 424         if (What.eq.2) call thermal(iprot)
 425 c        write (iout,*) "After thermal"
 426 c        call flush(iout)
 427       enddo ! iprot
 428 #ifdef MPI
 429       if (me1.eq.Master) then
 430 #endif
 431       do iprot=1,nprot
 432 #ifdef DEBUG
 433          write (iout,*) "iprot",iprot," elowest",
 434      &  ((elowest(ib,k,iprot),ib=1,nbeta(k,iprot)),k=1,natlike(iprot)+2)
 435         call flush(iout)
 436 #endif
 437         do ib=1,nbeta(2,iprot)
 438           fac=betaT(ib,2,iprot)
 439 c 6/15/05 AL Heat capacity
 440           zvtot(ib,iprot)=(esquare_ftot(ib,iprot)-
 441      &      emean_ftot(ib,iprot)**2)*fac*fac*Rgas
 442      &      -ebis_ftot(ib,iprot)+heatbase(iprot)
 443 #ifdef DEBUG
 444           write (iout,*) "iprot",iprot," ib",ib,
 445      &          " emeantot",emean_ftot(ib,iprot),
 446      &          " esquaretot",esquare_ftot(ib,iprot),
 447      &          " ebistot",ebis_ftot(ib,iprot),
 448      &          " zvtot",zvtot(ib,iprot)
 449 #endif
 450 #ifdef DEBUG
 451           write (iout,*) "beta",
 452      &        betaT(ib,2,iprot),
 453      &        " efree",efree(ib,2,iprot),
 454      &        " emean",emean_ftot(ib,iprot),
 455      &        " variance",esquare_ftot(ib,iprot)
 456           call flush(iout)
 457 #endif
 458 c 4/13/04 AL Correlation coefficients between energy and nativelikeness
 459 c      in the respective classes.
 460 #ifdef DEBUG
 461           write (iout,*) "iprot",iprot," ib",ib,
 462      &         " ii",ii
 463           write (iout,*) "emean",emean_ftot(ib,iprot),
 464      &         " esquare",esquare_ftot(ib,iprot)
 465 #endif
 466         enddo
 467       enddo
 468 #ifdef DEBUG
 469       write (iout,*)
 470      &   "Averages of nativelikeness and energy-nativelikeness",
 471      &          " correlation coefficients:"
 472       do iprot=1,nprot
 473       do i=1,natlike(iprot)
 474       do ib=1,nbeta(i+1,iprot)
 475         write (iout,*) "protein",iprot," property",i,
 476      &    " temperature",ib," natdim",natdim(i,iprot)," value(s)",
 477      &    (nuave(j,ib,i,iprot),j=1,natdim(i,iprot))
 478       enddo
 479       enddo
 480       enddo
 481       call flush(iout)
 482 #endif
 483       lprn=.false.
 484 #ifdef MPI
 485       endif
 486 #endif
 487 #ifdef DEBUG
 488         write (iout,*) "Calling CUTOFF_VIOLATION"
 489         call flush(iout)
 490 #endif
 491 c      write (iout,*) "CUTOFFEVAL",cutoffeval
 492       cutoffviol=.false.
 493       if (cutoffeval) call cutoff_violation
 494       t_func = t_func + tcpu() - tcpu_ini
 495 #ifdef DEBUG
 496         write (iout,*) "Exitting FUNC1"
 497         call flush(iout)
 498 #endif
 499       return
 500       end
 501 c-------------------------------------------------------------------------------
 502       subroutine targetfun(n,x,nf,f,uiparm,ufparm)
 503       implicit none
 504       include "DIMENSIONS"
 505       include "DIMENSIONS.ZSCOPT"
 506 #ifdef MPI
 507       include "mpif.h"
 508       integer IERROR,Errcode,status(MPI_STATUS_SIZE)
 509       include "COMMON.MPI"
 510 #endif
 511       integer n,nf
 512       double precision x(n),g(max_paropt)
 513       include "COMMON.WEIGHTS"
 514       include "COMMON.IOUNITS"
 515       include "COMMON.DERIV"
 516       include "COMMON.VMCPAR"
 517       include "COMMON.CLASSES"
 518       include "COMMON.GEO"
 519       include "COMMON.ENERGIES"
 520       include "COMMON.WEIGHTDER"
 521       include "COMMON.OPTIM"
 522       include "COMMON.TIME1"
 523       include "COMMON.COMPAR"
 524       include "COMMON.TORSION"
 525       double precision ff
 526       common /sru/ ff
 527       integer i,ii,it,inat,j,k,l,kk,iprot,ib
 528       logical bs
 529       external ufparm
 530       double precision urparm,rdif,maxlik_pdb
 531       integer uiparm
 532       double precision f,kara,viol,gnmr,gnmr1,aux
 533       logical lprn
 534       lprn = .true.
 535 #ifdef DEBUG
 536       write (iout,*) me,me1,"Calling TARGETFUN"
 537       write (iout,*) "n",n," x",(x(i),i=1,n)
 538 #endif
 539 c#undef DEBUG
 540       call write_restart(n,x)
 541       nfl=nf
 542       call x2w(n,x)
 543       f=0.0d0
 544       call func1(n,1,x)
 545 #ifdef NEWCORR
 546 c      if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 547       if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 548      &  fpmf = rdif(n,x,g,1)
 549       if (pdbpmf) fpdb = maxlik_pdb(n,x,g,1)
 550 c      endif
 551 #endif
 552 c Maximum likelihood contribution
 553       do iprot=1,nprot
 554         do iT=1,nbeta(1,iprot)
 555           f = f + weilik(iT,iprot)*sumlik(iT,iprot)
 556         enddo
 557       enddo
 558 #ifdef DEBUG
 559       write (iout,*)"target_fun: sumlik"
 560       do iprot=1,nprot
 561         write (iout,*) (sumlik(ii,iprot),ii=1,nbeta(1,iprot))
 562       enddo
 563       write (iout,*) "f before Cv =",f
 564       call flush(iout)
 565 #endif
 566       do iprot=1,nprot
 567 c 6/15/05 AL Contribution from Cv
 568         do ib=1,nbeta(2,iprot)
 569           f=f+weiCv(ib,iprot)*(zvtot(ib,iprot)-target_cv(ib,iprot))**2
 570         enddo
 571 c        write (iout,*) "target_fun from CV: f=",f
 572       enddo
 573 c      write (iout,*) "target_fun: f=",f
 574 c 8/2/2013 Contribution from conformation-dependent averages
 575       do iprot=1,nprot
 576         do inat=1,natlike(iprot)
 577           do ib=1,nbeta(inat+2,iprot)
 578             do k=1,natdim(inat,iprot)
 579               f = f+weinat(k,ib,inat,iprot)*
 580      &            (nuave(k,ib,inat,iprot)-nuexp(k,ib,inat,iprot))**2
 581 #ifdef DEBUG
 582               write (iout,*) "iprot",iprot," inat",inat,
 583      &        " ib",ib," k=",k," nuave",nuave(k,ib,inat,iprot),
 584      &        " nuexp",nuexp(k,ib,inat,iprot),
 585      &        " weight",weinat(k,ib,inat,iprot)
 586 #endif
 587             enddo ! k
 588           enddo ! ib
 589         enddo ! inat
 590       enddo ! iprot
 591 c      write (iout,*) "target_fun after adding nativelikeness: f=",f
 592       if (mode.lt.4) then
 593 c add constraints on weights
 594         kara=0.0d0
 595         do i=1,n
 596           kara=kara+gnmr1(x(i),x_low(i),x_up(i))
 597         enddo
 598         f=f+1.0d16*kara
 599       endif
 600       f = f + wpmf*fpmf + wpdbpmf*fpdb
 601       ff = f
 602 #ifdef DEBUG
 603       write(iout,*)"target_fun: f=",f,wpmf*fpmf,wpdbpmf*fpdb,1.0d16*kara
 604       call flush(iout)
 605 #endif
 606 c#undef DEBUG
 607       return
 608       end
 609 c-----------------------------------------------------------------------------
 610       subroutine targetgrad(n,x,nf,g,uiparm,rdum,ufparm)
 611       implicit none
 612       include "DIMENSIONS"
 613       include "DIMENSIONS.ZSCOPT"
 614       include "COMMON.TIME1"
 615 #ifdef MPI
 616       include "mpif.h"
 617       integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
 618 #endif
 619       integer n,nf
 620       double precision x(n),
 621      &  g(n),gg(max_paropt),ggg(max_paropt)
 622       include "COMMON.WEIGHTS"
 623       include "COMMON.ENERGIES"
 624       include "COMMON.CLASSES"
 625       include "COMMON.IOUNITS"
 626       include "COMMON.DERIV"
 627       include "COMMON.COMPAR"
 628       include "COMMON.WEIGHTDER"
 629       include "COMMON.OPTIM"
 630       include "COMMON.VMCPAR"
 631       include "COMMON.GEO"
 632       include "COMMON.NAMES"
 633       include "COMMON.VAR"
 634       include "COMMON.TORSION"
 635       double precision ff
 636       common /sru/ ff
 637 #ifdef MPI
 638       include "COMMON.MPI"
 639 #endif
 640       double precision xi,delta /1.0d-6/
 641       integer uiparm
 642       double precision urparm,rdif,maxlik_pdb
 643       external ufparm
 644       double precision zewn(MaxT,maxprot),aux,zscder
 645       double precision zewnvtot(maxT,maxprot)
 646       double precision sumlik0(maxT,maxprot),zv0tot(maxT,maxprot),
 647      &                 nuave0(maxdimnat,maxT,maxnatlike,maxprot)
 648       double precision geps(nntyp),gtor(maxtor_var),gang(maxang_var)
 649       double precision rdum
 650       double precision fac
 651       integer ll1,ll2,nl1,nl2,nn
 652       integer i,ii,ind,inat,j,k,kk,iprot,ib,n_weight
 653       integer iti,itj,jj,icant
 654       integer l,ll
 655       logical lprn
 656       double precision gnmrprim,gnmr1prim
 657       double precision tcpu_ini,tcpu
 658       logical in_sc_pair_list
 659       external in_sc_pair_list
 660       integer ind_shield /25/
 661       lprn=.false.
 662       tcpu_ini = tcpu()
 663       n_grad = n_grad+1
 664 #ifdef MPI
 665 #ifdef DEBUG
 666       write (iout,*) "Processor",me,me1," calling TARGETGRAD"
 667       call flush(iout)
 668 #endif
 669 #else
 670 #ifdef DEBUG
 671       write (iout,*) "Calling TARGETGRAD"
 672       call flush(iout)
 673 #endif
 674 #endif
 675 #ifdef DEBUG
 676       write (iout,*) "x",(x(i),i=1,n)
 677       call flush(iout)
 678 #endif
 679       icg=mod(nf,2)+1
 680       call x2w(n,x)
 681       call func1(n,1,x)
 682 #ifdef NEWCORR
 683 c AL 2/10/2018 Add PMF gradient
 684 c      if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 685       if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 686      &  aux = rdif(n,x,gg,2)
 687       if (pdbpmf) aux = maxlik_pdb(n,x,ggg,2)
 688 c        write (iout,*) "TARGETGRAD: After rdif"
 689 c        call flush(iout)
 690 c      endif
 691 #endif
 692       do iprot=1,nprot
 693         do ib=1,nbeta(2,iprot)
 694 c 6/15/05 AL Contributions from Cv
 695           zewnvtot(ib,iprot)=2*weiCv(ib,iprot)*
 696      &        (zvtot(ib,iprot)-target_cv(ib,iprot))
 697         enddo
 698       enddo
 699
 700       do k=1,n
 701         g(k)=0.0d0
 702       enddo
 703       if (mod_side) then
 704         do k=1,nntyp
 705           geps(k)=0.0d0
 706         enddo
 707       endif
 708       if (mod_tor .or. mod_sccor) then
 709         do k=1,maxtor_var
 710           gtor(k)=0.0d0
 711         enddo
 712       endif
 713       if (mod_ang) then
 714         do k=1,maxang_var
 715           gang(k)=0.0d0
 716         enddo
 717       endif
 718 c Maximum likelihood gradient
 719       do iprot=1,nprot
 720         do ib=1,nbeta(1,iprot)
 721         kk=0
 722         do k=1,n_ene
 723           if (imask(k).gt.0 .and. k.ne.ind_shield) then
 724             kk=kk+1
 725             g(kk)=g(kk)+weilik(ib,iprot)*sumlikder(kk,ib,iprot)
 726           endif
 727         enddo
 728         do k=1,nntyp
 729           geps(k)=geps(k)+weilik(ib,iprot)*sumlikeps(k,ib,iprot)
 730         enddo
 731         do k=1,ntor_var
 732           gtor(k)=gtor(k)+weilik(ib,iprot)*sumliktor(k,ib,iprot)
 733         enddo
 734         do k=1,nang_var
 735           gang(k)=gang(k)+weilik(ib,iprot)*sumlikang(k,ib,iprot)
 736         enddo
 737         enddo
 738       enddo
 739 #ifdef DEBUG
 740       write (iout,*) "gtor"
 741       do k=1,ntor_var
 742         write(iout,*) k,gtor(k)
 743       enddo
 744       write (iout,*) "nang_var",nang_var," gang"
 745       do k=1,nang_var
 746         write(iout,*) k,gang(k)
 747       enddo
 748 #endif
 749       call zderiv
 750 c Heat-capacity gradient
 751       do iprot=1,nprot
 752         do ib=1,nbeta(2,iprot)
 753           kk=0
 754           do k=1,n_ene
 755             if (imask(k).gt.0 .and. k.ne.ind_shield) then
 756               kk=kk+1
 757               g(kk)=g(kk)+zvdertot(kk,ib,iprot)*zewnvtot(ib,iprot)
 758             endif
 759           enddo
 760           do k=1,nntyp
 761             geps(k)=geps(k)+zewnvtot(ib,iprot)*zvepstot(k,ib,iprot)
 762           enddo
 763           do k=1,ntor_var
 764             gtor(k)=gtor(k)+zewnvtot(ib,iprot)*zvtortot(k,ib,iprot)
 765           enddo
 766           do k=1,nang_var
 767             gang(k)=gang(k)+zewnvtot(ib,iprot)*zvangtot(k,ib,iprot)
 768           enddo
 769         enddo
 770       enddo
 771 #ifdef DEBUG
 772       write (iout,*) "gtor"
 773       do k=1,ntor_var
 774         write(iout,*) k,gtor(k)
 775       enddo
 776       write (iout,*) "gang"
 777       do k=1,nang_var
 778         write(iout,*) k,gang(k)
 779       enddo
 780       write (iout,*) "grad: g",(g(i),i=1,n)
 781       call flush(iout)
 782 #endif
 783 c Derivatives of conformational-dependent properties
 784       do iprot=1,nprot
 785         do inat=1,natlike(iprot)
 786           do ib=1,nbeta(inat+2,iprot)
 787             call propderiv(ib,inat,iprot)
 788             do l=1,natdim(inat,iprot)
 789               aux=2*weinat(l,ib,inat,iprot)*
 790      &         (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
 791 c#ifdef DEBUG
 792 c              write (iout,*) "iprot",iprot," ib",ib,
 793 c     &         " inat",inat," l",l,
 794 c     &         " nuave",nuave(l,ib,inat,iprot),
 795 c     &         " weight",weinat(ib,inat,iprot),
 796 c     &         " zewnR",aux
 797 c              call flush(iout)
 798 c#endif
 799               kk=0
 800               do k=1,n_ene
 801                 if (imask(k).gt.0 .and. k.ne.ind_shield) then
 802                   kk=kk+1
 803                   g(kk)=g(kk)+aux*rder(kk,l)
 804                 endif
 805               enddo ! k
 806               do k=1,nntyp
 807                 geps(k)=geps(k)+aux*reps(k,l)
 808               enddo ! k
 809               do k=1,ntor_var
 810                 gtor(k)=gtor(k)+aux*rtor(k,l)
 811               enddo ! k
 812               do k=1,nang_var
 813                 gang(k)=gang(k)+aux*rang(k,l)
 814               enddo ! k
 815             enddo ! l
 816           enddo ! ib
 817         enddo
 818       enddo ! iprot
 819 #ifdef DEBUG
 820       write (iout,*) "gtor"
 821       do k=1,ntor_var
 822         write(iout,*) k,gtor(k)
 823       enddo
 824       write (iout,*) "gang"
 825       do k=1,nang_var
 826         write(iout,*) k,gang(k)
 827       enddo
 828 #endif
 829 c Compute numerically the remaining part of the gradient
 830 c      n_weight=0
 831 c      do i=1,n_ene
 832 c        if (imask(i).gt.0 .and. i.ne.ind_shield) n_weight=n_weight+1
 833 c      enddo
 834 c      ii=n_weight+17
 835       ii=n_weight_PMF
 836 c Sum up the gradient in SC epsilons, if needed
 837       if (mod_side) then
 838         do i=1,nsingle_sc
 839           ii=ii+1
 840           iti=ityp_ssc(i)
 841           g(ii)=g(ii)+geps(icant(iti,iti))
 842           do j=1,iti-1
 843             if (.not. in_sc_pair_list(iti,j)) then
 844               g(ii)=g(ii)+0.5d0*geps(icant(iti,j))
 845             endif
 846           enddo
 847         enddo
 848         do i=1,npair_sc
 849           ii=ii+1
 850           iti=ityp_psc(1,i)
 851           itj=ityp_psc(2,i)
 852           g(ii)=geps(icant(iti,itj))
 853         enddo
 854       endif
 855 c      write (iout,*) "n_ene",n_ene
 856       if (mod_tor .or. mod_sccor) then
 857         do i=1,ntor_var
 858           ii=ii+1
 859           g(ii)=gtor(i)
 860         enddo
 861       endif
 862       if (mod_ang) then
 863         do i=1,nang_var
 864           ii=ii+1
 865           g(ii)=gang(i)
 866         enddo
 867       endif
 868       n_weight=ii
 869 #ifdef DEBUG
 870       write (iout,*) "n",n," n_weight",n_weight
 871       write (iout,*) "After SC: grad:"
 872       do i=1,n
 873         write (iout,*) i,g(i)
 874       enddo
 875       call flush(iout)
 876 #endif
 877       call func1(n,1,x)
 878 #ifdef NEWCORR
 879 c AL 2/10/2018 Add PMF gradient
 880 c      if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 881       if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 882      &  aux = rdif(n,x,gg,0)
 883       if (pdbpmf) aux = maxlik_pdb(n,x,ggg,0)
 884 c        write (iout,*) "After rdif"
 885 c        call flush(iout)
 886 c      endif
 887 #endif
 888       do iprot=1,nprot
 889         do ib=1,nbeta(1,iprot)
 890           sumlik0(ib,iprot)=sumlik(ib,iprot)
 891         enddo
 892       enddo
 893       do iprot=1,nprot
 894         do ib=1,nbeta(2,iprot)
 895           zv0tot(ib,iprot)=zvtot(ib,iprot)
 896         enddo
 897       enddo
 898       do iprot=1,nprot
 899         do i=1,natlike(iprot)
 900           do ib=1,nbeta(i+2,iprot)
 901             do k=1,natdim(i,iprot)
 902               nuave0(k,ib,i,iprot)=nuave(k,ib,i,iprot)
 903             enddo
 904           enddo
 905         enddo
 906       enddo
 907       nn=n
 908       do iprot=1,nprot
 909         if (nbeta(2,iprot).gt.0) nn=nn-1
 910       enddo
 911       ii=nn
 912       do iprot=1,nprot
 913         if (nbeta(2,iprot).gt.0) then
 914           ii=ii+1
 915           do ib=1,nbeta(2,iprot)
 916             g(ii)=g(ii)+zewnvtot(ib,iprot)
 917           enddo
 918         endif
 919       enddo
 920 #ifdef DEBUG
 921       write (iout,*) "n",n," n_weight",n_weight
 922       write (iout,*) "After heat grad:"
 923       do i=1,n
 924         write (iout,*) i,g(i)
 925       enddo
 926       call flush(iout)
 927 #endif
 928 c      write (iout,*) "n_weight",n_weight," n",n
 929       do i=n_weight+1,nn
 930         xi=x(i)
 931         x(i)=xi+delta
 932 c        write (iout,*) "i",i," x",x(i),xi
 933         call x2w(n,x)
 934         call func1(n,1,x)
 935 #ifdef NEWCORR
 936 c AL 2/10/2018 Add PMF gradient
 937 c      if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 938       if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 939      &  aux = rdif(n,x,gg,0)
 940       if (pdbpmf) aux = maxlik_pdb(n,x,ggg,0)
 941 c        write (iout,*) "After rdif"
 942 c        call flush(iout)
 943 c      endif
 944 #endif
 945 c Maximum likelihood
 946         do iprot=1,nprot
 947           do ib=1,nbeta(1,iprot)
 948             g(i)=g(i)+weilik(ib,iprot)*(sumlik(ib,iprot)
 949      &             -sumlik0(ib,iprot))/delta
 950 c            write (iout,*) "iprot",iprot," g",g(i)
 951           enddo
 952         enddo
 953 #ifdef DEBUG
 954         write (iout,*) "init numerical derivatives"
 955         call flush(iout)
 956 #endif
 957 c Contribution from Cv
 958         do iprot=1,nprot
 959           do ib=1,nbeta(2,iprot)
 960             g(i)=g(i)+zewnvtot(ib,iprot)*
 961      &      (zvtot(ib,iprot)-zv0tot(ib,iprot))/delta
 962           enddo
 963         enddo
 964 c Contribution from average nativelikenss and correlation coefficients
 965 c between energy and nativelikeness
 966         do iprot=1,nprot
 967           do ii=1,natlike(iprot)
 968             do ib=1,nbeta(ii+2,iprot)
 969               do k=1,natdim(i,iprot)
 970                aux=2*weinat(k,ib,inat,iprot)*
 971      &          (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
 972                g(i)=g(i)+aux*(nuave(k,ib,ii,iprot)-
 973      &           nuave0(k,ib,ii,iprot))/delta
 974               enddo ! k
 975             enddo
 976           enddo
 977         enddo
 978         x(i)=xi
 979       enddo
 980 #ifdef DEBUG
 981       write (iout,*) "n",n," n_weight",n_weight
 982       write (iout,*) "After num grad:"
 983       call flush(iout)
 984       do k=1,n
 985         write (iout,*) k,g(k)
 986       enddo
 987       call flush(iout)
 988 #endif
 989       call func1(n,1,x)
 990 #ifdef NEWCORR
 991 c AL 2/10/2018 Add PMF gradient
 992 c      if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 993       if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 994      &  aux = rdif(n,x,gg,2)
 995       if (pdbpmf)   aux = maxlik_pdb(n,x,ggg,0)
 996 c        write (iout,*) "After rdif"
 997 c        call flush(iout)
 998 c      endif
 999 #endif
1000 #ifdef DEBUG
1001         write (iout,*) "finish numerical derivatives"
1002         call flush(iout)
1003 #endif
1004 c
1005 c 1/10/2007 AL Correction: ENDIF...IF insterted to prevent unmatched
1006 c           SEND is I_LOCAL_CHECK.GT.0
1007 c
1008 #ifdef DEBUG
1009       write (iout,*) "mode",mode
1010 #endif
1011 #ifdef NEWCORR
1012 c AL 2/10/2018 Add PMF gradient and PDB-PMF-maxlik gradinet
1013       if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
1014         do i=1,n
1015           g(i)=g(i)-wpmf*gg(i)+wpdbpmf*ggg(i)
1016         enddo
1017       endif
1018 #endif
1019       if (mode.gt.3) then
1020 c transform the gradient to the x space
1021         do i=1,n
1022           g(i)=g(i)*(x_up(i)-x_low(i))/(pi*(1.0d0+(x(i)
1023      &     -0.5d0*(x_up(i)+x_low(i)))**2))
1024         enddo
1025       else
1026 c add constraints on weights
1027 c        write (iout,*) "Gradient in constraints"
1028         do i=1,n
1029           g(i)=g(i)+1.0d16*gnmr1prim(x(i),x_low(i),x_up(i))
1030 #ifdef DEBUG
1031           write (iout,*) i,x(i),x_low(i),x_up(i),
1032      &      gnmr1prim(x(i),x_low(i),x_up(i))
1033 #endif
1034         enddo
1035       endif
1036 #ifdef DEBUG
1037         write (iout,*) "grad"
1038         do i=1,n
1039           write (iout,*) i,g(i)
1040         enddo
1041 #endif
1042       t_grad = t_grad + tcpu() - tcpu_ini
1043 #ifdef DEBUG
1044       write (iout,*) "Exit targetgrad"
1045       call flush(iout)
1046 #endif
1047       return
1048       end
1049 c------------------------------------------------------------------------------
1050       subroutine matout(ndim,n,a)
1051       implicit real*8 (a-h,o-z)
1052       include "DIMENSIONS.ZSCOPT"
1053       include "COMMON.IOUNITS"
1054       double precision a(ndim)
1055       character*6 line6 / '------' /
1056       character*12 line12 / '------------' /
1057       dimension b(8)
1058       do 1 i=1,n,6
1059       nlim=min0(i+5,n)
1060       write (iout,1000) (k,k=i,nlim)
1061       write (iout,1020) line6,(line12,k=i,nlim)
1062  1000 format (/7x,6(i7,5x))
1063  1020 format (a6,6a12)
1064       do 2 j=i,n
1065       jlim=min0(j,nlim)
1066       jlim1=jlim-i+1
1067       do 3 k=i,jlim
1068     3 b(k-i+1)=a(icant(i,j))
1069       write (iout,1010) j,(b(k),k=1,jlim1)
1070     2 continue
1071     1 continue
1072  1010 format (i3,3x,6(1pe11.4,1x))
1073       return
1074       end
1075 c------------------------------------------------------------------------------
1076       subroutine zderiv
1077       implicit none
1078       include "DIMENSIONS"
1079       include "DIMENSIONS.ZSCOPT"
1080 #ifdef MPI
1081       include "mpif.h"
1082       integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1083 #endif
1084       integer n,nf
1085       include "COMMON.WEIGHTS"
1086       include "COMMON.ENERGIES"
1087       include "COMMON.IOUNITS"
1088       include "COMMON.DERIV"
1089       include "COMMON.WEIGHTDER"
1090       include "COMMON.VMCPAR"
1091       include "COMMON.GEO"
1092       include "COMMON.NAMES"
1093       include "COMMON.VAR"
1094       include "COMMON.CLASSES"
1095       include "COMMON.THERMAL"
1096 #ifdef MPI
1097       include "COMMON.MPI"
1098 #endif
1099       include "COMMON.TIME1"
1100       double precision aux,zscder
1101       double precision fac,fac2
1102       double precision efj(2),emj(2),varj(2),edecum,edecum1,ebisdecum,
1103      &  e2decum
1104       double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1105      & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1106      & enepsmixsqfj(nntyp,2),eavpfjprim(max_ene,2),
1107      & entoravefj(0:3,maxtor_var,2),entoravefjprim(0:3,maxtor_var,2),
1108      & entormixfj(0:3,maxtor_var,2),entormixsqfj(0:3,maxtor_var,2),
1109      & enangavefj(0:3,maxang_var,2),enangavefjprim(0:3,maxang_var,2),
1110      & enangmixfj(0:3,maxang_var,2),enangmixsqfj(0:3,maxang_var,2)
1111       integer i,ii,j,k,kk,iprot,ib
1112       integer iti,itj,jj,icant
1113       integer l,ll
1114       logical lprn
1115       integer ind_shield /25/
1116       lprn=.false.
1117       do iprot=1,nprot
1118 C Heat capacity
1119         do ib=1,nbeta(2,iprot)
1120           fac = betaT(ib,2,iprot)
1121           fac2 = fac*fac
1122           kk=0
1123           do k=1,n_ene
1124             if (imask(k).gt.0 .and. k.ne.ind_shield) then
1125               kk=kk+1
1126               ebisdecum=emix_pfbistot(k,ib,iprot)-
1127      &         eave_pftot(k,ib,iprot)
1128      &         *ebis_ftot(ib,iprot)
1129               edecum=emix_pfprimtot(k,ib,iprot)
1130      &         -eave_pfprimtot(k,ib,iprot)
1131      &         *emean_ftot(ib,iprot)
1132               edecum1=emix_pftot(k,ib,iprot)
1133      &         -eave_pftot(k,ib,iprot)
1134      &         *emean_ftot(ib,iprot)
1135               e2decum=emixsq_pftot(k,ib,iprot)
1136      &         -eave_pftot(k,ib,iprot)
1137      &         *esquare_ftot(ib,iprot)
1138               zvdertot(kk,ib,iprot)=Rgas*fac2*(2*edecum
1139      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1140      &         *edecum1))
1141      &         -eave_pfbistot(k,ib,iprot)
1142      &         +fac*ebisdecum
1143 #ifdef DEBUG
1144               write (iout,*) "k=",k,
1145      &            "eave_pftot",eave_pftot(k,ib,iprot),
1146      &            " emix_pftot",emix_pftot(k,ib,iprot),
1147      &            " emix_pfprimtot",emix_pfprimtot(k,ib,iprot),
1148      &            " emix_pfbistot",emix_pfbistot(k,ib,iprot),
1149      &            " eave_pfprimtot",eave_pfprimtot(k,ib,iprot),
1150      &            " eave_pfbistot",eave_pfbistot(k,ib,iprot),
1151      &            " emixsq_pftot",emixsq_pftot(k,ib,iprot),
1152      &            " ebis_ftot",ebis_ftot(ib,iprot),
1153      &            " emean_ftot",emean_ftot(ib,iprot)
1154               write (iout,*) "k=",k," edecum",edecum,
1155      &            " edecum1",edecum1,
1156      &            " e2decum",e2decum," ebisdecum",ebisdecum,
1157      &            " zvdertot",
1158      &             zvdertot(kk,ib,iprot)
1159 #endif
1160             endif
1161           enddo
1162           if (mod_side) then
1163             do k=1,nntyp
1164               ebisdecum=eneps_mix_fbistot(k,ib,iprot)-
1165      &          enepsave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1166               edecum=eneps_mix_ftot(k,ib,iprot)
1167      &         -enepsave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1168               e2decum=eneps_mixsq_ftot(k,ib,iprot)
1169      &         -enepsave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1170               zvepstot(k,ib,iprot)=ww(1)*(Rgas*fac2*(2*edecum
1171      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1172      &         *edecum))+fac*ebisdecum)
1173 #ifdef DEBUG
1174               write (iout,*) ib,k,"ebisdecum",ebisdecum,
1175      &         " edecum",edecum," e2decum",e2decum," zvepstot",
1176      &          zvepstot(k,ibb,iprot)
1177               write (iout,*) k,
1178      &         " enepsave_ftot",enepsave_ftot(k,ib,iprot),
1179      &         " eneps_mix_ftot",eneps_mix_ftot(k,ib,iprot),
1180      &         " eneps_mix_fbistot",eneps_mix_fbistot(k,ib,iprot),
1181      &         " eneps_mixsq_ftot",eneps_mixsq_ftot(k,ib,iprot),
1182      &         " emean_ftot",emean_ftot(ib,iprot),
1183      &         " ebis_ftot",ebis_ftot(ib,iprot)
1184 #endif
1185             enddo
1186           endif
1187           if (mod_tor .or. mod_sccor) then
1188             do k=1,ntor_var
1189               ebisdecum=entor_mix_fbistot(k,ib,iprot)-
1190      &          entorave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1191               edecum=entor_mix_fprimtot(k,ib,iprot)
1192      &         -entorave_fprimtot(k,ib,iprot)*emean_ftot(ib,iprot)
1193               edecum1=entor_mix_ftot(k,ib,iprot)
1194      &         -entorave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1195               e2decum=entor_mixsq_ftot(k,ib,iprot)
1196      &         -entorave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1197               zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1198      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1199      &         *edecum1))-entorave_fbistot(k,ib,iprot)
1200      &         +fac*ebisdecum
1201 #ifdef DEBUG
1202               write (iout,*) k,
1203      &         " entorave_ftot",entorave_ftot(k,ib,iprot),
1204      &         " entor_mix_ftot",entor_mix_ftot(k,ib,iprot),
1205      &         " entor_mix_fprimtot",entor_mix_fprimtot(k,ib,iprot),
1206      &         " entor_mix_fbistot",entor_mix_fbistot(k,ib,iprot),
1207      &         " entorave_fprimtot",entorave_fprimtot(k,ib,iprot),
1208      &         " entorave_fbistot",entorave_fbistot(k,ib,iprot),
1209      &         " entor_mixsq_ftot",entor_mixsq_ftot(k,ib,iprot),
1210      &         " emean_ftot",emean_ftot(ib,iprot),
1211      &         " ebis_ftot",ebis_ftot(ib,iprot)
1212               write (iout,*) ib,k," ebisdecum",ebisdecum,
1213      &         " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1214      &         " zvtortot",zvtortot(k,ibb,iprot)
1215 #endif
1216             enddo
1217           endif
1218           if (mod_ang) then
1219             do k=1,nang_var
1220               ebisdecum=0.0d0
1221               edecum=0.0d0
1222               edecum1=enang_mix_ftot(k,ib,iprot)
1223      &         -enangave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1224               e2decum=enang_mixsq_ftot(k,ib,iprot)
1225      &         -enangave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1226               zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1227      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1228      &         *edecum1))-entorave_fbistot(k,ib,iprot)
1229      &         +fac*ebisdecum
1230 #ifdef DEBUG
1231               write (iout,*) k,
1232      &         " enangave_ftot",enangave_ftot(k,ib,iprot),
1233      &         " enang_mix_ftot",enang_mix_ftot(k,ib,iprot),
1234      &         " enang_mixsq_ftot",enang_mixsq_ftot(k,ib,iprot),
1235      &         " emean_ftot",emean_ftot(ib,iprot)
1236               write (iout,*) ib,k," ebisdecum",ebisdecum,
1237      &         " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1238      &         " zvtortot",zvangtot(k,ibb,iprot)
1239 #endif
1240             enddo
1241           endif
1242 #ifdef DEBUG
1243           write (iout,*) "iprot",iprot,
1244      &        " ib",ib
1245           write (iout,*) "zvder",
1246      &       (zvdertot(k,ibb,iprot),k=1,kk)
1247 #endif
1248         enddo ! ib
1249       enddo ! iprot
1250       return
1251       end
1252 c------------------------------------------------------------------------------
1253       subroutine propderiv(ib,inat,iprot)
1254 c Compute the derivatives of conformation-dependent averages in force-field
1255 c parameters.
1256       implicit none
1257       include "DIMENSIONS"
1258       include "DIMENSIONS.ZSCOPT"
1259 #ifdef MPI
1260       include "mpif.h"
1261       integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1262 #endif
1263       integer n,nf,inat
1264       include "COMMON.WEIGHTS"
1265       include "COMMON.ENERGIES"
1266       include "COMMON.CLASSES"
1267       include "COMMON.IOUNITS"
1268       include "COMMON.DERIV"
1269       include "COMMON.WEIGHTDER"
1270       include "COMMON.VMCPAR"
1271       include "COMMON.GEO"
1272       include "COMMON.OPTIM"
1273       include "COMMON.NAMES"
1274       include "COMMON.COMPAR"
1275 #ifdef MPI
1276       include "COMMON.MPI"
1277 #endif
1278       include "COMMON.TIME1"
1279       include "COMMON.VAR"
1280       double precision aux,zscder
1281       double precision fac
1282       double precision efj(2),emj(2),varj(2)
1283       double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1284      & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1285      & enepsmixsqfj(nntyp,2)
1286       integer ll1,nl1
1287       integer i,ii,ind,j,k,kk,iprot,ib,ibb,l1
1288       integer iti,itj,jj,icant
1289       integer l,ll
1290       logical lprn
1291       integer ind_shield /25/
1292 c      write (iout,*) "Entered PROPDERIV"
1293       lprn=.false.
1294       fac = betaT(ib,inat+2,iprot)
1295 c        write (iout,*) "derivatives: iprot",iprot,
1296 c     &     " ib",ib," l",l
1297       do l=1,natdim(inat,iprot)
1298         kk=0
1299         do k=1,n_ene
1300           if (imask(k).gt.0 .and. k.ne.ind_shield) then
1301             kk=kk+1
1302             rder(kk,l)=-fac*(nu_pf(k,l,ib,inat,iprot)
1303      &       -nuave(l,ib,inat,iprot)*
1304      &       eave_nat_pftot(k,ib,inat,iprot))
1305           endif
1306         enddo
1307       enddo
1308       if (mod_side) then
1309         do l=1,natdim(inat,iprot)
1310           do k=1,nntyp
1311             reps(k,l)=-ww(1)*fac*(nuepsave_f(k,l,ib,inat,iprot)-
1312      &        nuave(l,ib,inat,iprot)*enepsave_nat_ftot(k,ib,inat,iprot))
1313           enddo
1314         enddo
1315       endif
1316       if (mod_tor .or. mod_sccor) then
1317         do l=1,natdim(inat,iprot)
1318           do k=1,ntor_var
1319             rtor(k,l)=-fac*(nutorave_f(k,l,ib,inat,iprot)-
1320      &        nuave(l,inat,ib,iprot)*entorave_nat_ftot(k,ib,inat,iprot))
1321           enddo
1322         enddo
1323       endif
1324       if (mod_ang) then
1325         do l=1,natdim(inat,iprot)
1326           do k=1,nang_var
1327             rang(k,l)=-fac*(nuangave_f(k,l,ib,inat,iprot)-
1328      &        nuave(l,inat,ib,iprot)*enangave_nat_ftot(k,ib,inat,iprot))
1329           enddo
1330         enddo
1331       endif
1332       return
1333       end
1334 c-------------------------------------------------------------------------------
1335       subroutine ene_eval(iprot,i,ii)
1336       implicit none
1337       include "DIMENSIONS"
1338       include "DIMENSIONS.ZSCOPT"
1339 #ifdef MPI
1340       include "mpif.h"
1341       integer IERROR,ErrCode
1342       include "COMMON.MPI"
1343 #endif
1344       include "COMMON.CONTROL"
1345       include "COMMON.COMPAR"
1346       include "COMMON.WEIGHTS"
1347       include "COMMON.WEIGHTDER"
1348       include "COMMON.CLASSES"
1349       include "COMMON.ENERGIES"
1350       include "COMMON.IOUNITS"
1351       include "COMMON.VMCPAR"
1352       include "COMMON.NAMES"
1353       include "COMMON.INTERACT"
1354       include "COMMON.TIME1"
1355       include "COMMON.CHAIN"
1356       include "COMMON.VAR"
1357       include "COMMON.GEO"
1358       include "COMMON.SCCOR"
1359       include "COMMON.TORSION"
1360       include "COMMON.LOCAL"
1361       include "COMMON.FFIELD"
1362 C Define local variables
1363       integer i,ii,ij,jj,j,k,ik,jk,l,ll,lll,iprot,iblock
1364       double precision energia(0:max_ene)
1365       double precision etoti,enepsjk
1366       double precision ftune_eps
1367       external ftune_eps
1368       logical lprn
1369       integer icant
1370       external icant
1371       integer ind_shield /25/
1372       lprn=.true.
1373       ii = ii+1
1374 #ifdef DEBUG
1375       write (iout,*) "ene_eval: iprot",iprot," i",i," ii",ii
1376 #endif
1377       if (init_ene .or. (mod_fourier(nloctyp).gt.0
1378      &  .or. mod_fouriertor(nloctyp).gt.0
1379      &  .or. mod_elec .or. mod_scp .or. mod_tor .or. mod_sccor
1380      &  .or. mod_ang .or. imask(ind_shield).gt.0) ) then
1381 #ifdef DEBUG
1382         write (iout,*) "FUNC1: doing UNRES"
1383         write (iout,*) "Protein",iprot," i",i," ii",ii," nres",nres
1384         call flush(iout)
1385 #endif
1386         call restore_ccoords(iprot,ii)
1387         call int_from_cart1(.false.)
1388 #ifdef DEBUG
1389         write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1390         write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1391         write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1392         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1393         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1394         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1395         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1396         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1397
1398 c        call intout
1399 c        call cartprint
1400 #endif
1401         call etotal(energia(0))
1402 c        if (natlike(iprot).gt.0)
1403 c     &    call natsimil(i,ii,iprot,.false.)
1404         e_total(i,iprot)=energia(0)
1405         do j=1,n_ene
1406           enetb(ii,j,iprot)=energia(j)
1407         enddo
1408 c#define DEBUG
1409 #ifdef DEBUG
1410         write (iout,*) "ww",ww(:n_ene)
1411         write (iout,'(2i5,18(1pe12.4))') i,ii,
1412      &  (enetb(ii,j,iprot),j=1,n_ene)
1413         write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1414      &       i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)
1415         call enerprint(energia(0))
1416 #endif
1417 c#undef DEBUG
1418         if (mod_side) then
1419           write (iout,*) "FUN nntyp",nntyp," eneps",eneps_temp(1,45)
1420           do j=1,nntyp
1421             do k=1,2
1422               eneps(k,j,ii,iprot)=eneps_temp(k,j)
1423             enddo
1424           enddo
1425 c#ifdef DEBUG
1426 c            write (iout,*) "Matrix eneps"
1427 c            do j=1,ntyp
1428 c              do k=1,j
1429 c                write (iout,*) j,k,icant(j,k),
1430 c     &           eneps_temp(1,icant(j,k)),
1431 c     &           eneps_temp(2,icant(j,k))
1432 c              enddo
1433 c            enddo
1434 c            call flush(iout)
1435 c#endif
1436         endif
1437         if (mod_tor .or. mod_sccor) then
1438 c          write (iout,*) "etor",enetb(ii,13,iprot),enetb(ii,19,iprot)
1439           k=0
1440           do iblock=1,2
1441           do ij=-ntortyp+1,ntortyp-1
1442             do jj=-ntortyp+1,ntortyp-1
1443               if (mask_tor(0,jj,ij,iblock).eq.1) then
1444                 k=k+1
1445                 entor(k,ii,iprot)=etor_temp(0,0,jj,ij,iblock)
1446 c                write (iout,*) "etor_temp",restyp(itortyp(j)),
1447 c     &            restyp(itortyp(i)),
1448 c     &            k,etor_temp(ji,ii)
1449               else if (mask_tor(0,jj,ij,iblock).gt.1) then
1450                 if (tor_mode.eq.0) then
1451                 do l=1,nterm(jj,ij,iblock)
1452                   k=k+1
1453                   entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1454                 enddo
1455                 if (mask_tor(0,jj,ij,iblock).gt.2) then
1456                   do l=nterm(jj,ij,iblock)+1,2*nterm(jj,ij,iblock)
1457                     k=k+1
1458                     entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1459                   enddo
1460                 endif
1461                 else if (tor_mode.eq.1) then
1462
1463                 do l=1,nterm_kcc(jj,ij)
1464                   do ll=1,nterm_kcc_Tb(jj,ij)
1465                     do lll=1,nterm_kcc_Tb(jj,ij)
1466                       k=k+1
1467                       entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1468                     enddo ! lll
1469                   enddo ! ll
1470                 enddo
1471                 if (mask_tor(0,jj,ij,1).gt.2) then
1472                   do l=nterm_kcc(jj,ij)+1,2*nterm_kcc(jj,ij)
1473                     do ll=1,nterm_kcc_Tb(jj,ij)
1474                       do lll=1,nterm_kcc_Tb(jj,ij)
1475                         k=k+1
1476                         entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1477                       enddo
1478                     enddo
1479                   enddo
1480
1481                 endif
1482                 endif
1483               endif
1484             enddo
1485           enddo
1486           enddo
1487           do ij=-nsccortyp+1,nsccortyp-1
1488             do jj=-nsccortyp+1,nsccortyp-1
1489               do ll=1,3
1490               if (mask_tor(ll,jj,ij,1).eq.1) then
1491                 k=k+1
1492                 entor(k,ii,iprot)=etor_temp(0,ll,jj,ij,1)
1493 c                write (iout,*) "etor_temp",restyp(isccortyp(jj)),
1494 c     &            restyp(isccortyp(ij)),
1495 c     &            k,etor_temp(jj,ij)
1496               else if (mask_tor(ll,jj,ij,1).gt.1) then
1497                 do l=1,nterm_sccor(jj,ij)
1498                   k=k+1
1499                   entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1500                 enddo
1501                 do l=nterm_sccor(jj,ij)+1,2*nterm_sccor(jj,ij)
1502                   k=k+1
1503                   entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1504                 enddo
1505               endif
1506               enddo
1507             enddo
1508           enddo
1509         endif
1510 c AL 8/7/17: optimization of the bending potentials
1511         k=0
1512         if (mod_ang) then
1513           do ij=0,nthetyp-1
1514             if (mask_ang(ij).eq.0) cycle
1515             do jj=0,nbend_kcc_Tb(ij)
1516               k=k+1
1517               enbend(k,ii,iprot)=ebend_temp_kcc(jj,ij)*wang
1518 c              write (iout,*) "ij",ij," jj",jj," k",k," ebend_temp_kcc",
1519 c     &         ebend_temp_kcc(jj,ij)
1520             enddo
1521           enddo
1522         endif
1523 c#define DEBUG
1524 #ifdef DEBUG
1525         write (iout,*) "Conformation",i
1526 c        write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1527 c        write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1528 c        write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1529 c        write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1530 c        write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1531 c        write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1532 c        write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1533 c        write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1534         call enerprint(energia(0))
1535 c        write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1536 c     &       i,energia(0),eini(i,iprot),entfac(i,iprot)
1537 c        write (iout,'(i5,18(1pe12.4))') ii,
1538 c     &  (enetb(ii,j,iprot),j=1,n_ene)
1539 c        write (iout,'(i5,16(1pe12.4))') i,
1540 c     &  (energia(j),j=1,n_ene),energia(0)
1541 c        print *,"Processor",me,me1," computing energies",i
1542 #endif
1543 c#undef DEBUG
1544         if (energia(0).ge.1.0d99) then
1545           write (iout,*) "FUNC1:CHUJ NASTAPIL in energy evaluation for",
1546      &  " point",i,". Probably NaNs in some of the energy components."
1547           write (iout,*) "The components of the energy are:"
1548           call enerprint(energia(0))
1549           write (iout,*) "Conformation",i,ii
1550           write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1551           write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1552           write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1553           write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1554           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1555           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1556           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1557           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1558           write (iout,*) "Calculation terminated at this point.",
1559      &       " Check the database of conformations"
1560           call flush(iout)
1561 #ifdef MPI
1562           call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
1563 #endif
1564           stop "SEVERE error in energy calculation"
1565         endif
1566       else
1567         if (mod_side) then
1568           enetb(ii,1,iprot)=0.0d0
1569           do j=1,ntyp
1570             do k=1,j
1571               enetb(ii,1,iprot)=enetb(ii,1,iprot)+
1572      &            ftune_eps(eps(j,k))*eneps(1,icant(j,k),ii,iprot)+
1573      &            eps(j,k)*eneps(2,icant(j,k),ii,iprot)
1574             enddo
1575           enddo
1576         endif
1577         etoti=0.0d0
1578         do j=1,n_ene
1579           etoti=etoti+enetb(ii,j,iprot)*ww(j)
1580         enddo
1581 c#define DEBUG
1582 #ifdef DEBUG
1583         write (iout,*) "i",i," ww",ww
1584         write (iout,*) "enetb",(enetb(ii,j,iprot),j=1,n_ene)
1585 #endif
1586 c#undef DEBUG
1587         e_total(i,iprot)=etoti
1588 c        write (iout,*) "FUNC1 natlike",iprot,natlike(iprot)
1589         if (natlike(iprot).gt.0) then
1590           call restore_ccoords(iprot,ii)
1591           call int_from_cart1(.false.)
1592 c          call natsimil(i,ii,iprot,.false.)
1593 c          write (iout,*)
1594 c     &   "natsimil",i,ii,(nu(k,ii,iprot),k=1,natlike(iprot))
1595         endif
1596 #ifdef DEBUG
1597         write (iout,'(i5,18(1pe12.4))') i,
1598      &  (enetb(ii,j,iprot),j=1,n_ene)
1599         write (iout,'(18(1pe12.4))')
1600      &  (eneps(1,j,ii,iprot),j=201,210)
1601         write(iout,'(4hENER,i10,3f15.3)')
1602      &       i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)
1603 #endif
1604       endif
1605       return
1606       end