source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/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 #ifdef DEBUG
 597       write (iout,*) "target_fun: f=",f,wpmf*fpmf,1.0d16*kara
 598       call flush(iout)
 599 #endif
 600 c#undef DEBUG
 601       return
 602       end
 603 c-----------------------------------------------------------------------------
 604       subroutine targetgrad(n,x,nf,g,uiparm,rdum,ufparm)
 605       implicit none
 606       include "DIMENSIONS"
 607       include "DIMENSIONS.ZSCOPT"
 608       include "COMMON.TIME1"
 609 #ifdef MPI
 610       include "mpif.h"
 611       integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
 612 #endif
 613       integer n,nf
 614       double precision x(n),
 615      &  g(n),gg(max_paropt)
 616       include "COMMON.WEIGHTS"
 617       include "COMMON.ENERGIES"
 618       include "COMMON.CLASSES"
 619       include "COMMON.IOUNITS"
 620       include "COMMON.DERIV"
 621       include "COMMON.COMPAR"
 622       include "COMMON.WEIGHTDER"
 623       include "COMMON.OPTIM"
 624       include "COMMON.VMCPAR"
 625       include "COMMON.GEO"
 626       include "COMMON.NAMES"
 627       include "COMMON.VAR"
 628       include "COMMON.TORSION"
 629       double precision ff
 630       common /sru/ ff
 631 #ifdef MPI
 632       include "COMMON.MPI"
 633 #endif
 634       double precision xi,delta /1.0d-6/
 635       integer uiparm
 636       double precision urparm,rdif
 637       external ufparm
 638       double precision zewn(MaxT,maxprot),aux,zscder
 639       double precision zewnvtot(maxT,maxprot)
 640       double precision sumlik0(maxT,maxprot),zv0tot(maxT,maxprot),
 641      &                 nuave0(maxdimnat,maxT,maxnatlike,maxprot)
 642       double precision geps(nntyp),gtor(maxtor_var),gang(maxang_var)
 643       double precision rdum
 644       double precision fac
 645       integer ll1,ll2,nl1,nl2,nn
 646       integer i,ii,ind,inat,j,k,kk,iprot,ib,n_weight
 647       integer iti,itj,jj,icant
 648       integer l,ll
 649       logical lprn
 650       double precision gnmrprim,gnmr1prim
 651       double precision tcpu_ini,tcpu
 652       logical in_sc_pair_list
 653       external in_sc_pair_list
 654       integer ind_shield /25/
 655       lprn=.false.
 656       tcpu_ini = tcpu()
 657       n_grad = n_grad+1
 658 #ifdef MPI
 659 #ifdef DEBUG
 660       write (iout,*) "Processor",me,me1," calling TARGETGRAD"
 661       call flush(iout)
 662 #endif
 663 #else
 664 #ifdef DEBUG
 665       write (iout,*) "Calling TARGETGRAD"
 666       call flush(iout)
 667 #endif
 668 #endif
 669 #ifdef DEBUG
 670       write (iout,*) "x",(x(i),i=1,n)
 671       call flush(iout)
 672 #endif
 673       icg=mod(nf,2)+1
 674       call x2w(n,x)
 675       call func1(n,1,x)
 676 #ifdef NEWCORR
 677 c AL 2/10/2018 Add PMF gradient
 678       if (mod_fourier(nloctyp).gt.0) then
 679         aux = rdif(n,x,gg,2)
 680 c        write (iout,*) "After rdif"
 681 c        call flush(iout)
 682       endif
 683 #endif
 684       do iprot=1,nprot
 685         do ib=1,nbeta(2,iprot)
 686 c 6/15/05 AL Contributions from Cv
 687           zewnvtot(ib,iprot)=2*weiCv(ib,iprot)*
 688      &        (zvtot(ib,iprot)-target_cv(ib,iprot))
 689         enddo
 690       enddo
 691
 692       do k=1,n
 693         g(k)=0.0d0
 694       enddo
 695       if (mod_side) then
 696         do k=1,nntyp
 697           geps(k)=0.0d0
 698         enddo
 699       endif
 700       if (mod_tor .or. mod_sccor) then
 701         do k=1,maxtor_var
 702           gtor(k)=0.0d0
 703         enddo
 704       endif
 705       if (mod_ang) then
 706         do k=1,maxang_var
 707           gang(k)=0.0d0
 708         enddo
 709       endif
 710 c Maximum likelihood gradient
 711       do iprot=1,nprot
 712         do ib=1,nbeta(1,iprot)
 713         kk=0
 714         do k=1,n_ene
 715           if (imask(k).gt.0 .and. k.ne.ind_shield) then
 716             kk=kk+1
 717             g(kk)=g(kk)+weilik(ib,iprot)*sumlikder(kk,ib,iprot)
 718           endif
 719         enddo
 720         do k=1,nntyp
 721           geps(k)=geps(k)+weilik(ib,iprot)*sumlikeps(k,ib,iprot)
 722         enddo
 723         do k=1,ntor_var
 724           gtor(k)=gtor(k)+weilik(ib,iprot)*sumliktor(k,ib,iprot)
 725         enddo
 726         do k=1,nang_var
 727           gang(k)=gang(k)+weilik(ib,iprot)*sumlikang(k,ib,iprot)
 728         enddo
 729         enddo
 730       enddo
 731 #ifdef DEBUG
 732       write (iout,*) "gtor"
 733       do k=1,ntor_var
 734         write(iout,*) k,gtor(k)
 735       enddo
 736       write (iout,*) "nang_var",nang_var," gang"
 737       do k=1,nang_var
 738         write(iout,*) k,gang(k)
 739       enddo
 740 #endif
 741       call zderiv
 742 c Heat-capacity gradient
 743       do iprot=1,nprot
 744         do ib=1,nbeta(2,iprot)
 745           kk=0
 746           do k=1,n_ene
 747             if (imask(k).gt.0 .and. k.ne.ind_shield) then
 748               kk=kk+1
 749               g(kk)=g(kk)+zvdertot(kk,ib,iprot)*zewnvtot(ib,iprot)
 750             endif
 751           enddo
 752           do k=1,nntyp
 753             geps(k)=geps(k)+zewnvtot(ib,iprot)*zvepstot(k,ib,iprot)
 754           enddo
 755           do k=1,ntor_var
 756             gtor(k)=gtor(k)+zewnvtot(ib,iprot)*zvtortot(k,ib,iprot)
 757           enddo
 758           do k=1,nang_var
 759             gang(k)=gang(k)+zewnvtot(ib,iprot)*zvangtot(k,ib,iprot)
 760           enddo
 761         enddo
 762       enddo
 763 #ifdef DEBUG
 764       write (iout,*) "gtor"
 765       do k=1,ntor_var
 766         write(iout,*) k,gtor(k)
 767       enddo
 768       write (iout,*) "gang"
 769       do k=1,nang_var
 770         write(iout,*) k,gang(k)
 771       enddo
 772       write (iout,*) "grad: g",(g(i),i=1,n)
 773       call flush(iout)
 774 #endif
 775 c Derivatives of conformational-dependent properties
 776       do iprot=1,nprot
 777         do inat=1,natlike(iprot)
 778           do ib=1,nbeta(inat+2,iprot)
 779             call propderiv(ib,inat,iprot)
 780             do l=1,natdim(inat,iprot)
 781               aux=2*weinat(l,ib,inat,iprot)*
 782      &         (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
 783 #ifdef DEBUG
 784               write (iout,*) "iprot",iprot," ib",ib,
 785      &         " inat",inat," l",l,
 786      &         " nuave",nuave(l,ib,inat,iprot),
 787      &         " weight",weinat(ib,inat,iprot),
 788      &         " zewnR",aux
 789               call flush(iout)
 790 #endif
 791               kk=0
 792               do k=1,n_ene
 793                 if (imask(k).gt.0 .and. k.ne.ind_shield) then
 794                   kk=kk+1
 795                   g(kk)=g(kk)+aux*rder(kk,l)
 796                 endif
 797               enddo ! k
 798               do k=1,nntyp
 799                 geps(k)=geps(k)+aux*reps(k,l)
 800               enddo ! k
 801               do k=1,ntor_var
 802                 gtor(k)=gtor(k)+aux*rtor(k,l)
 803               enddo ! k
 804               do k=1,nang_var
 805                 gang(k)=gang(k)+aux*rang(k,l)
 806               enddo ! k
 807             enddo ! l
 808           enddo ! ib
 809         enddo
 810       enddo ! iprot
 811 #ifdef DEBUG
 812       write (iout,*) "gtor"
 813       do k=1,ntor_var
 814         write(iout,*) k,gtor(k)
 815       enddo
 816       write (iout,*) "gang"
 817       do k=1,nang_var
 818         write(iout,*) k,gang(k)
 819       enddo
 820 #endif
 821 c Compute numerically the remaining part of the gradient
 822 c      n_weight=0
 823 c      do i=1,n_ene
 824 c        if (imask(i).gt.0 .and. i.ne.ind_shield) n_weight=n_weight+1
 825 c      enddo
 826 c      ii=n_weight+17
 827       ii=n_weight_PMF
 828 c Sum up the gradient in SC epsilons, if needed
 829       if (mod_side) then
 830         do i=1,nsingle_sc
 831           ii=ii+1
 832           iti=ityp_ssc(i)
 833           g(ii)=g(ii)+geps(icant(iti,iti))
 834           do j=1,iti-1
 835             if (.not. in_sc_pair_list(iti,j)) then
 836               g(ii)=g(ii)+0.5d0*geps(icant(iti,j))
 837             endif
 838           enddo
 839         enddo
 840         do i=1,npair_sc
 841           ii=ii+1
 842           iti=ityp_psc(1,i)
 843           itj=ityp_psc(2,i)
 844           g(ii)=geps(icant(iti,itj))
 845         enddo
 846       endif
 847 c      write (iout,*) "n_ene",n_ene
 848       if (mod_tor .or. mod_sccor) then
 849         do i=1,ntor_var
 850           ii=ii+1
 851           g(ii)=gtor(i)
 852         enddo
 853       endif
 854       if (mod_ang) then
 855         do i=1,nang_var
 856           ii=ii+1
 857           g(ii)=gang(i)
 858         enddo
 859       endif
 860       n_weight=ii
 861 #ifdef DEBUG
 862       write (iout,*) "n",n," n_weight",n_weight
 863       write (iout,*) "After SC: grad:"
 864       do i=1,n
 865         write (iout,*) i,g(i)
 866       enddo
 867       call flush(iout)
 868 #endif
 869       call func1(n,1,x)
 870 #ifdef NEWCORR
 871 c AL 2/10/2018 Add PMF gradient
 872       if (mod_fourier(nloctyp).gt.0) then
 873         aux = rdif(n,x,gg,0)
 874 c        write (iout,*) "After rdif"
 875 c        call flush(iout)
 876       endif
 877 #endif
 878       do iprot=1,nprot
 879         do ib=1,nbeta(1,iprot)
 880           sumlik0(ib,iprot)=sumlik(ib,iprot)
 881         enddo
 882       enddo
 883       do iprot=1,nprot
 884         do ib=1,nbeta(2,iprot)
 885           zv0tot(ib,iprot)=zvtot(ib,iprot)
 886         enddo
 887       enddo
 888       do iprot=1,nprot
 889         do i=1,natlike(iprot)
 890           do ib=1,nbeta(i+2,iprot)
 891             do k=1,natdim(i,iprot)
 892               nuave0(k,ib,i,iprot)=nuave(k,ib,i,iprot)
 893             enddo
 894           enddo
 895         enddo
 896       enddo
 897       nn=n
 898       do iprot=1,nprot
 899         if (nbeta(2,iprot).gt.0) nn=nn-1
 900       enddo
 901       ii=nn
 902       do iprot=1,nprot
 903         if (nbeta(2,iprot).gt.0) then
 904           ii=ii+1
 905           do ib=1,nbeta(2,iprot)
 906             g(ii)=g(ii)+zewnvtot(ib,iprot)
 907           enddo
 908         endif
 909       enddo
 910 #ifdef DEBUG
 911       write (iout,*) "n",n," n_weight",n_weight
 912       write (iout,*) "After heat grad:"
 913       do i=1,n
 914         write (iout,*) i,g(i)
 915       enddo
 916       call flush(iout)
 917 #endif
 918 c      write (iout,*) "n_weight",n_weight," n",n
 919       do i=n_weight+1,nn
 920         xi=x(i)
 921         x(i)=xi+delta
 922 c        write (iout,*) "i",i," x",x(i),xi
 923         call x2w(n,x)
 924         call func1(n,1,x)
 925 #ifdef NEWCORR
 926 c AL 2/10/2018 Add PMF gradient
 927       if (mod_fourier(nloctyp).gt.0) then
 928         aux = rdif(n,x,gg,0)
 929 c        write (iout,*) "After rdif"
 930 c        call flush(iout)
 931       endif
 932 #endif
 933 c Maximum likelihood
 934         do iprot=1,nprot
 935           do ib=1,nbeta(1,iprot)
 936             g(i)=g(i)+weilik(ib,iprot)*(sumlik(ib,iprot)
 937      &             -sumlik0(ib,iprot))/delta
 938 c            write (iout,*) "iprot",iprot," g",g(i)
 939           enddo
 940         enddo
 941 #ifdef DEBUG
 942         write (iout,*) "init numerical derivatives"
 943         call flush(iout)
 944 #endif
 945 c Contribution from Cv
 946         do iprot=1,nprot
 947           do ib=1,nbeta(2,iprot)
 948             g(i)=g(i)+zewnvtot(ib,iprot)*
 949      &      (zvtot(ib,iprot)-zv0tot(ib,iprot))/delta
 950           enddo
 951         enddo
 952 c Contribution from average nativelikenss and correlation coefficients
 953 c between energy and nativelikeness
 954         do iprot=1,nprot
 955           do ii=1,natlike(iprot)
 956             do ib=1,nbeta(ii+2,iprot)
 957               do k=1,natdim(i,iprot)
 958                aux=2*weinat(k,ib,inat,iprot)*
 959      &          (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
 960                g(i)=g(i)+aux*(nuave(k,ib,ii,iprot)-
 961      &           nuave0(k,ib,ii,iprot))/delta
 962               enddo ! k
 963             enddo
 964           enddo
 965         enddo
 966         x(i)=xi
 967       enddo
 968 #ifdef DEBUG
 969       write (iout,*) "n",n," n_weight",n_weight
 970       write (iout,*) "After num grad:"
 971       do k=1,n
 972         write (iout,*) k,g(k)
 973       enddo
 974       call flush(iout)
 975 #endif
 976       call func1(n,1,x)
 977 #ifdef NEWCORR
 978 c AL 2/10/2018 Add PMF gradient
 979       if (mod_fourier(nloctyp).gt.0) then
 980         aux = rdif(n,x,gg,2)
 981 c        write (iout,*) "After rdif"
 982 c        call flush(iout)
 983       endif
 984 #endif
 985 #ifdef DEBUG
 986         write (iout,*) "finish numerical derivatives"
 987         call flush(iout)
 988 #endif
 989 c
 990 c 1/10/2007 AL Correction: ENDIF...IF insterted to prevent unmatched
 991 c           SEND is I_LOCAL_CHECK.GT.0
 992 c
 993 #ifdef DEBUG
 994       write (iout,*) "mode",mode
 995 #endif
 996 #ifdef NEWCORR
 997 c AL 2/10/2018 Add PMF gradient
 998       if (mod_fourier(nloctyp).gt.0) then
 999         do i=1,n
1000           g(i)=g(i)-wpmf*gg(i)
1001         enddo
1002       endif
1003 #endif
1004       if (mode.gt.3) then
1005 c transform the gradient to the x space
1006         do i=1,n
1007           g(i)=g(i)*(x_up(i)-x_low(i))/(pi*(1.0d0+(x(i)
1008      &     -0.5d0*(x_up(i)+x_low(i)))**2))
1009         enddo
1010       else
1011 c add constraints on weights
1012 c        write (iout,*) "Gradient in constraints"
1013         do i=1,n
1014           g(i)=g(i)+1.0d16*gnmr1prim(x(i),x_low(i),x_up(i))
1015 #ifdef DEBUG
1016           write (iout,*) i,x(i),x_low(i),x_up(i),
1017      &      gnmr1prim(x(i),x_low(i),x_up(i))
1018 #endif
1019         enddo
1020       endif
1021 #ifdef DEBUG
1022         write (iout,*) "grad"
1023         do i=1,n
1024           write (iout,*) i,g(i)
1025         enddo
1026 #endif
1027       t_grad = t_grad + tcpu() - tcpu_ini
1028 #ifdef DEBUG
1029       write (iout,*) "Exit targetgrad"
1030       call flush(iout)
1031 #endif
1032       return
1033       end
1034 c------------------------------------------------------------------------------
1035       subroutine matout(ndim,n,a)
1036       implicit real*8 (a-h,o-z)
1037       include "DIMENSIONS.ZSCOPT"
1038       include "COMMON.IOUNITS"
1039       double precision a(ndim)
1040       character*6 line6 / '------' /
1041       character*12 line12 / '------------' /
1042       dimension b(8)
1043       do 1 i=1,n,6
1044       nlim=min0(i+5,n)
1045       write (iout,1000) (k,k=i,nlim)
1046       write (iout,1020) line6,(line12,k=i,nlim)
1047  1000 format (/7x,6(i7,5x))
1048  1020 format (a6,6a12)
1049       do 2 j=i,n
1050       jlim=min0(j,nlim)
1051       jlim1=jlim-i+1
1052       do 3 k=i,jlim
1053     3 b(k-i+1)=a(icant(i,j))
1054       write (iout,1010) j,(b(k),k=1,jlim1)
1055     2 continue
1056     1 continue
1057  1010 format (i3,3x,6(1pe11.4,1x))
1058       return
1059       end
1060 c------------------------------------------------------------------------------
1061       subroutine zderiv
1062       implicit none
1063       include "DIMENSIONS"
1064       include "DIMENSIONS.ZSCOPT"
1065 #ifdef MPI
1066       include "mpif.h"
1067       integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1068 #endif
1069       integer n,nf
1070       include "COMMON.WEIGHTS"
1071       include "COMMON.ENERGIES"
1072       include "COMMON.IOUNITS"
1073       include "COMMON.DERIV"
1074       include "COMMON.WEIGHTDER"
1075       include "COMMON.VMCPAR"
1076       include "COMMON.GEO"
1077       include "COMMON.NAMES"
1078       include "COMMON.VAR"
1079       include "COMMON.CLASSES"
1080       include "COMMON.THERMAL"
1081 #ifdef MPI
1082       include "COMMON.MPI"
1083 #endif
1084       include "COMMON.TIME1"
1085       double precision aux,zscder
1086       double precision fac,fac2
1087       double precision efj(2),emj(2),varj(2),edecum,edecum1,ebisdecum,
1088      &  e2decum
1089       double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1090      & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1091      & enepsmixsqfj(nntyp,2),eavpfjprim(max_ene,2),
1092      & entoravefj(0:3,maxtor_var,2),entoravefjprim(0:3,maxtor_var,2),
1093      & entormixfj(0:3,maxtor_var,2),entormixsqfj(0:3,maxtor_var,2),
1094      & enangavefj(0:3,maxang_var,2),enangavefjprim(0:3,maxang_var,2),
1095      & enangmixfj(0:3,maxang_var,2),enangmixsqfj(0:3,maxang_var,2)
1096       integer i,ii,j,k,kk,iprot,ib
1097       integer iti,itj,jj,icant
1098       integer l,ll
1099       logical lprn
1100       integer ind_shield /25/
1101       lprn=.false.
1102       do iprot=1,nprot
1103 C Heat capacity
1104         do ib=1,nbeta(2,iprot)
1105           fac = betaT(ib,2,iprot)
1106           fac2 = fac*fac
1107           kk=0
1108           do k=1,n_ene
1109             if (imask(k).gt.0 .and. k.ne.ind_shield) then
1110               kk=kk+1
1111               ebisdecum=emix_pfbistot(k,ib,iprot)-
1112      &         eave_pftot(k,ib,iprot)
1113      &         *ebis_ftot(ib,iprot)
1114               edecum=emix_pfprimtot(k,ib,iprot)
1115      &         -eave_pfprimtot(k,ib,iprot)
1116      &         *emean_ftot(ib,iprot)
1117               edecum1=emix_pftot(k,ib,iprot)
1118      &         -eave_pftot(k,ib,iprot)
1119      &         *emean_ftot(ib,iprot)
1120               e2decum=emixsq_pftot(k,ib,iprot)
1121      &         -eave_pftot(k,ib,iprot)
1122      &         *esquare_ftot(ib,iprot)
1123               zvdertot(kk,ib,iprot)=Rgas*fac2*(2*edecum
1124      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1125      &         *edecum1))
1126      &         -eave_pfbistot(k,ib,iprot)
1127      &         +fac*ebisdecum
1128 #ifdef DEBUG
1129               write (iout,*) "k=",k,
1130      &            "eave_pftot",eave_pftot(k,ib,iprot),
1131      &            " emix_pftot",emix_pftot(k,ib,iprot),
1132      &            " emix_pfprimtot",emix_pfprimtot(k,ib,iprot),
1133      &            " emix_pfbistot",emix_pfbistot(k,ib,iprot),
1134      &            " eave_pfprimtot",eave_pfprimtot(k,ib,iprot),
1135      &            " eave_pfbistot",eave_pfbistot(k,ib,iprot),
1136      &            " emixsq_pftot",emixsq_pftot(k,ib,iprot),
1137      &            " ebis_ftot",ebis_ftot(ib,iprot),
1138      &            " emean_ftot",emean_ftot(ib,iprot)
1139               write (iout,*) "k=",k," edecum",edecum,
1140      &            " edecum1",edecum1,
1141      &            " e2decum",e2decum," ebisdecum",ebisdecum,
1142      &            " zvdertot",
1143      &             zvdertot(kk,ib,iprot)
1144 #endif
1145             endif
1146           enddo
1147           if (mod_side) then
1148             do k=1,nntyp
1149               ebisdecum=eneps_mix_fbistot(k,ib,iprot)-
1150      &          enepsave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1151               edecum=eneps_mix_ftot(k,ib,iprot)
1152      &         -enepsave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1153               e2decum=eneps_mixsq_ftot(k,ib,iprot)
1154      &         -enepsave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1155               zvepstot(k,ib,iprot)=ww(1)*(Rgas*fac2*(2*edecum
1156      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1157      &         *edecum))+fac*ebisdecum)
1158 #ifdef DEBUG
1159               write (iout,*) ib,k,"ebisdecum",ebisdecum,
1160      &         " edecum",edecum," e2decum",e2decum," zvepstot",
1161      &          zvepstot(k,ibb,iprot)
1162               write (iout,*) k,
1163      &         " enepsave_ftot",enepsave_ftot(k,ib,iprot),
1164      &         " eneps_mix_ftot",eneps_mix_ftot(k,ib,iprot),
1165      &         " eneps_mix_fbistot",eneps_mix_fbistot(k,ib,iprot),
1166      &         " eneps_mixsq_ftot",eneps_mixsq_ftot(k,ib,iprot),
1167      &         " emean_ftot",emean_ftot(ib,iprot),
1168      &         " ebis_ftot",ebis_ftot(ib,iprot)
1169 #endif
1170             enddo
1171           endif
1172           if (mod_tor .or. mod_sccor) then
1173             do k=1,ntor_var
1174               ebisdecum=entor_mix_fbistot(k,ib,iprot)-
1175      &          entorave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1176               edecum=entor_mix_fprimtot(k,ib,iprot)
1177      &         -entorave_fprimtot(k,ib,iprot)*emean_ftot(ib,iprot)
1178               edecum1=entor_mix_ftot(k,ib,iprot)
1179      &         -entorave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1180               e2decum=entor_mixsq_ftot(k,ib,iprot)
1181      &         -entorave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1182               zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1183      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1184      &         *edecum1))-entorave_fbistot(k,ib,iprot)
1185      &         +fac*ebisdecum
1186 #ifdef DEBUG
1187               write (iout,*) k,
1188      &         " entorave_ftot",entorave_ftot(k,ib,iprot),
1189      &         " entor_mix_ftot",entor_mix_ftot(k,ib,iprot),
1190      &         " entor_mix_fprimtot",entor_mix_fprimtot(k,ib,iprot),
1191      &         " entor_mix_fbistot",entor_mix_fbistot(k,ib,iprot),
1192      &         " entorave_fprimtot",entorave_fprimtot(k,ib,iprot),
1193      &         " entorave_fbistot",entorave_fbistot(k,ib,iprot),
1194      &         " entor_mixsq_ftot",entor_mixsq_ftot(k,ib,iprot),
1195      &         " emean_ftot",emean_ftot(ib,iprot),
1196      &         " ebis_ftot",ebis_ftot(ib,iprot)
1197               write (iout,*) ib,k," ebisdecum",ebisdecum,
1198      &         " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1199      &         " zvtortot",zvtortot(k,ibb,iprot)
1200 #endif
1201             enddo
1202           endif
1203           if (mod_ang) then
1204             do k=1,nang_var
1205               ebisdecum=0.0d0
1206               edecum=0.0d0
1207               edecum1=enang_mix_ftot(k,ib,iprot)
1208      &         -enangave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1209               e2decum=enang_mixsq_ftot(k,ib,iprot)
1210      &         -enangave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1211               zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1212      &         -fac*(e2decum-2*emean_ftot(ib,iprot)
1213      &         *edecum1))-entorave_fbistot(k,ib,iprot)
1214      &         +fac*ebisdecum
1215 #ifdef DEBUG
1216               write (iout,*) k,
1217      &         " enangave_ftot",enangave_ftot(k,ib,iprot),
1218      &         " enang_mix_ftot",enang_mix_ftot(k,ib,iprot),
1219      &         " enang_mixsq_ftot",enang_mixsq_ftot(k,ib,iprot),
1220      &         " emean_ftot",emean_ftot(ib,iprot)
1221               write (iout,*) ib,k," ebisdecum",ebisdecum,
1222      &         " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1223      &         " zvtortot",zvangtot(k,ibb,iprot)
1224 #endif
1225             enddo
1226           endif
1227 #ifdef DEBUG
1228           write (iout,*) "iprot",iprot,
1229      &        " ib",ib
1230           write (iout,*) "zvder",
1231      &       (zvdertot(k,ibb,iprot),k=1,kk)
1232 #endif
1233         enddo ! ib
1234       enddo ! iprot
1235       return
1236       end
1237 c------------------------------------------------------------------------------
1238       subroutine propderiv(ib,inat,iprot)
1239 c Compute the derivatives of conformation-dependent averages in force-field
1240 c parameters.
1241       implicit none
1242       include "DIMENSIONS"
1243       include "DIMENSIONS.ZSCOPT"
1244 #ifdef MPI
1245       include "mpif.h"
1246       integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1247 #endif
1248       integer n,nf,inat
1249       include "COMMON.WEIGHTS"
1250       include "COMMON.ENERGIES"
1251       include "COMMON.CLASSES"
1252       include "COMMON.IOUNITS"
1253       include "COMMON.DERIV"
1254       include "COMMON.WEIGHTDER"
1255       include "COMMON.VMCPAR"
1256       include "COMMON.GEO"
1257       include "COMMON.OPTIM"
1258       include "COMMON.NAMES"
1259       include "COMMON.COMPAR"
1260 #ifdef MPI
1261       include "COMMON.MPI"
1262 #endif
1263       include "COMMON.TIME1"
1264       include "COMMON.VAR"
1265       double precision aux,zscder
1266       double precision fac
1267       double precision efj(2),emj(2),varj(2)
1268       double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1269      & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1270      & enepsmixsqfj(nntyp,2)
1271       integer ll1,nl1
1272       integer i,ii,ind,j,k,kk,iprot,ib,ibb,l1
1273       integer iti,itj,jj,icant
1274       integer l,ll
1275       logical lprn
1276       integer ind_shield /25/
1277 c      write (iout,*) "Entered PROPDERIV"
1278       lprn=.false.
1279       fac = betaT(ib,inat+2,iprot)
1280 c        write (iout,*) "derivatives: iprot",iprot,
1281 c     &     " ib",ib," l",l
1282       do l=1,natdim(inat,iprot)
1283         kk=0
1284         do k=1,n_ene
1285           if (imask(k).gt.0 .and. k.ne.ind_shield) then
1286             kk=kk+1
1287             rder(kk,l)=-fac*(nu_pf(k,l,ib,inat,iprot)
1288      &       -nuave(l,ib,inat,iprot)*
1289      &       eave_nat_pftot(k,ib,inat,iprot))
1290           endif
1291         enddo
1292       enddo
1293       if (mod_side) then
1294         do l=1,natdim(inat,iprot)
1295           do k=1,nntyp
1296             reps(k,l)=-ww(1)*fac*(nuepsave_f(k,l,ib,inat,iprot)-
1297      &        nuave(l,ib,inat,iprot)*enepsave_nat_ftot(k,ib,inat,iprot))
1298           enddo
1299         enddo
1300       endif
1301       if (mod_tor .or. mod_sccor) then
1302         do l=1,natdim(inat,iprot)
1303           do k=1,ntor_var
1304             rtor(k,l)=-fac*(nutorave_f(k,l,ib,inat,iprot)-
1305      &        nuave(l,inat,ib,iprot)*entorave_nat_ftot(k,ib,inat,iprot))
1306           enddo
1307         enddo
1308       endif
1309       if (mod_ang) then
1310         do l=1,natdim(inat,iprot)
1311           do k=1,nang_var
1312             rang(k,l)=-fac*(nuangave_f(k,l,ib,inat,iprot)-
1313      &        nuave(l,inat,ib,iprot)*enangave_nat_ftot(k,ib,inat,iprot))
1314           enddo
1315         enddo
1316       endif
1317       return
1318       end
1319 c-------------------------------------------------------------------------------
1320       subroutine ene_eval(iprot,i,ii)
1321       implicit none
1322       include "DIMENSIONS"
1323       include "DIMENSIONS.ZSCOPT"
1324 #ifdef MPI
1325       include "mpif.h"
1326       integer IERROR,ErrCode
1327       include "COMMON.MPI"
1328 #endif
1329       include "COMMON.CONTROL"
1330       include "COMMON.COMPAR"
1331       include "COMMON.WEIGHTS"
1332       include "COMMON.WEIGHTDER"
1333       include "COMMON.CLASSES"
1334       include "COMMON.ENERGIES"
1335       include "COMMON.IOUNITS"
1336       include "COMMON.VMCPAR"
1337       include "COMMON.NAMES"
1338       include "COMMON.INTERACT"
1339       include "COMMON.TIME1"
1340       include "COMMON.CHAIN"
1341       include "COMMON.VAR"
1342       include "COMMON.GEO"
1343       include "COMMON.SCCOR"
1344       include "COMMON.TORSION"
1345       include "COMMON.LOCAL"
1346       include "COMMON.FFIELD"
1347 C Define local variables
1348       integer i,ii,ij,jj,j,k,ik,jk,l,ll,lll,iprot,iblock
1349       double precision energia(0:max_ene)
1350       double precision etoti,enepsjk
1351       double precision ftune_eps
1352       external ftune_eps
1353       logical lprn
1354       integer icant
1355       external icant
1356       integer ind_shield /25/
1357       lprn=.true.
1358       ii = ii+1
1359 #ifdef DEBUG
1360       write (iout,*) "ene_eval: iprot",iprot," i",i," ii",ii
1361 #endif
1362       if (init_ene .or. (mod_fourier(nloctyp).gt.0
1363      &  .or. mod_elec .or. mod_scp .or. mod_tor .or. mod_sccor
1364      &  .or. mod_ang .or. imask(ind_shield).gt.0) ) then
1365 #ifdef DEBUG
1366         write (iout,*) "FUNC1: doing UNRES"
1367         write (iout,*) "Protein",iprot," i",i," ii",ii," nres",nres
1368         call flush(iout)
1369 #endif
1370         call restore_ccoords(iprot,ii)
1371         call int_from_cart1(.false.)
1372 #ifdef DEBUG
1373         write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1374         write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1375         write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1376         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1377         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1378         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1379         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1380         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1381
1382 c        call intout
1383 c        call cartprint
1384 #endif
1385         call etotal(energia(0))
1386 c        if (natlike(iprot).gt.0)
1387 c     &    call natsimil(i,ii,iprot,.false.)
1388         e_total(i,iprot)=energia(0)
1389         do j=1,n_ene
1390           enetb(ii,j,iprot)=energia(j)
1391         enddo
1392 c#define DEBUG
1393 #ifdef DEBUG
1394         write (iout,*) "ww",ww(:n_ene)
1395         write (iout,'(2i5,18(1pe12.4))') i,ii,
1396      &  (enetb(ii,j,iprot),j=1,n_ene)
1397         write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1398      &       i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)
1399         call enerprint(energia(0))
1400 #endif
1401 c#undef DEBUG
1402         if (mod_side) then
1403           write (iout,*) "FUN nntyp",nntyp," eneps",eneps_temp(1,45)
1404           do j=1,nntyp
1405             do k=1,2
1406               eneps(k,j,ii,iprot)=eneps_temp(k,j)
1407             enddo
1408           enddo
1409 c#ifdef DEBUG
1410 c            write (iout,*) "Matrix eneps"
1411 c            do j=1,ntyp
1412 c              do k=1,j
1413 c                write (iout,*) j,k,icant(j,k),
1414 c     &           eneps_temp(1,icant(j,k)),
1415 c     &           eneps_temp(2,icant(j,k))
1416 c              enddo
1417 c            enddo
1418 c            call flush(iout)
1419 c#endif
1420         endif
1421         if (mod_tor .or. mod_sccor) then
1422 c          write (iout,*) "etor",enetb(ii,13,iprot),enetb(ii,19,iprot)
1423           k=0
1424           do iblock=1,2
1425           do ij=-ntortyp+1,ntortyp-1
1426             do jj=-ntortyp+1,ntortyp-1
1427               if (mask_tor(0,jj,ij,iblock).eq.1) then
1428                 k=k+1
1429                 entor(k,ii,iprot)=etor_temp(0,0,jj,ij,iblock)
1430 c                write (iout,*) "etor_temp",restyp(itortyp(j)),
1431 c     &            restyp(itortyp(i)),
1432 c     &            k,etor_temp(ji,ii)
1433               else if (mask_tor(0,jj,ij,iblock).gt.1) then
1434                 if (tor_mode.eq.0) then
1435                 do l=1,nterm(jj,ij,iblock)
1436                   k=k+1
1437                   entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1438                 enddo
1439                 if (mask_tor(0,jj,ij,iblock).gt.2) then
1440                   do l=nterm(jj,ij,iblock)+1,2*nterm(jj,ij,iblock)
1441                     k=k+1
1442                     entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1443                   enddo
1444                 endif
1445                 else if (tor_mode.eq.1) then
1446
1447                 do l=1,nterm_kcc(jj,ij)
1448                   do ll=1,nterm_kcc_Tb(jj,ij)
1449                     do lll=1,nterm_kcc_Tb(jj,ij)
1450                       k=k+1
1451                       entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1452                     enddo ! lll
1453                   enddo ! ll
1454                 enddo
1455                 if (mask_tor(0,jj,ij,1).gt.2) then
1456                   do l=nterm_kcc(jj,ij)+1,2*nterm_kcc(jj,ij)
1457                     do ll=1,nterm_kcc_Tb(jj,ij)
1458                       do lll=1,nterm_kcc_Tb(jj,ij)
1459                         k=k+1
1460                         entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1461                       enddo
1462                     enddo
1463                   enddo
1464
1465                 endif
1466                 endif
1467               endif
1468             enddo
1469           enddo
1470           enddo
1471           do ij=-nsccortyp+1,nsccortyp-1
1472             do jj=-nsccortyp+1,nsccortyp-1
1473               do ll=1,3
1474               if (mask_tor(ll,jj,ij,1).eq.1) then
1475                 k=k+1
1476                 entor(k,ii,iprot)=etor_temp(0,ll,jj,ij,1)
1477 c                write (iout,*) "etor_temp",restyp(isccortyp(jj)),
1478 c     &            restyp(isccortyp(ij)),
1479 c     &            k,etor_temp(jj,ij)
1480               else if (mask_tor(ll,jj,ij,1).gt.1) then
1481                 do l=1,nterm_sccor(jj,ij)
1482                   k=k+1
1483                   entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1484                 enddo
1485                 do l=nterm_sccor(jj,ij)+1,2*nterm_sccor(jj,ij)
1486                   k=k+1
1487                   entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1488                 enddo
1489               endif
1490               enddo
1491             enddo
1492           enddo
1493         endif
1494 c AL 8/7/17: optimization of the bending potentials
1495         k=0
1496         if (mod_ang) then
1497           do ij=0,nthetyp-1
1498             if (mask_ang(ij).eq.0) cycle
1499             do jj=1,nbend_kcc_Tb(ij)
1500               k=k+1
1501               enbend(k,ii,iprot)=ebend_temp_kcc(jj,ij)*wang
1502 c              write (iout,*) "ij",ij," jj",jj," k",k," ebend_temp_kcc",
1503 c     &         ebend_temp_kcc(jj,ij)
1504             enddo
1505           enddo
1506         endif
1507 c#define DEBUG
1508 #ifdef DEBUG
1509         write (iout,*) "Conformation",i
1510 c        write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1511 c        write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1512 c        write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1513 c        write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1514 c        write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1515 c        write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1516 c        write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1517 c        write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1518         call enerprint(energia(0))
1519 c        write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1520 c     &       i,energia(0),eini(i,iprot),entfac(i,iprot)
1521 c        write (iout,'(i5,18(1pe12.4))') ii,
1522 c     &  (enetb(ii,j,iprot),j=1,n_ene)
1523 c        write (iout,'(i5,16(1pe12.4))') i,
1524 c     &  (energia(j),j=1,n_ene),energia(0)
1525 c        print *,"Processor",me,me1," computing energies",i
1526 #endif
1527 c#undef DEBUG
1528         if (energia(0).ge.1.0d99) then
1529           write (iout,*) "FUNC1:CHUJ NASTAPIL in energy evaluation for",
1530      &  " point",i,". Probably NaNs in some of the energy components."
1531           write (iout,*) "The components of the energy are:"
1532           call enerprint(energia(0))
1533           write (iout,*) "Conformation",i,ii
1534           write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1535           write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1536           write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1537           write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1538           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1539           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1540           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1541           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1542           write (iout,*) "Calculation terminated at this point.",
1543      &       " Check the database of conformations"
1544           call flush(iout)
1545 #ifdef MPI
1546           call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
1547 #endif
1548           stop "SEVERE error in energy calculation"
1549         endif
1550       else
1551         if (mod_side) then
1552           enetb(ii,1,iprot)=0.0d0
1553           do j=1,ntyp
1554             do k=1,j
1555               enetb(ii,1,iprot)=enetb(ii,1,iprot)+
1556      &            ftune_eps(eps(j,k))*eneps(1,icant(j,k),ii,iprot)+
1557      &            eps(j,k)*eneps(2,icant(j,k),ii,iprot)
1558             enddo
1559           enddo
1560         endif
1561         etoti=0.0d0
1562         do j=1,n_ene
1563           etoti=etoti+enetb(ii,j,iprot)*ww(j)
1564         enddo
1565 c#define DEBUG
1566 #ifdef DEBUG
1567         write (iout,*) "i",i," ww",ww
1568         write (iout,*) "enetb",(enetb(ii,j,iprot),j=1,n_ene)
1569 #endif
1570 c#undef DEBUG
1571         e_total(i,iprot)=etoti
1572 c        write (iout,*) "FUNC1 natlike",iprot,natlike(iprot)
1573         if (natlike(iprot).gt.0) then
1574           call restore_ccoords(iprot,ii)
1575           call int_from_cart1(.false.)
1576 c          call natsimil(i,ii,iprot,.false.)
1577 c          write (iout,*)
1578 c     &   "natsimil",i,ii,(nu(k,ii,iprot),k=1,natlike(iprot))
1579         endif
1580 #ifdef DEBUG
1581         write (iout,'(i5,18(1pe12.4))') i,
1582      &  (enetb(ii,j,iprot),j=1,n_ene)
1583         write (iout,'(18(1pe12.4))')
1584      &  (eneps(1,j,ii,iprot),j=201,210)
1585         write(iout,'(4hENER,i10,3f15.3)')
1586      &       i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)
1587 #endif
1588       endif
1589       return
1590       end