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