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