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