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