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

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