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