source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/minimize.F

   1       subroutine minimize_init(n,iv,v,d)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       integer liv,lv
   6       parameter (liv=60,lv=(77+max_paropt*(max_paropt+17)/2))
   7       include "COMMON.OPTIM"
   8       include "COMMON.MINPAR"
   9       include "COMMON.IOUNITS"
  10       include "COMMON.TIME1"
  11       integer n,iv(liv)
  12       double precision d(max_paropt),v(1:lv)
  13       integer i
  14
  15       llocal=.false.
  16       call deflt(2,iv,liv,lv,v)
  17 * 12 means fresh start, dont call deflt
  18       iv(1)=12
  19 * max num of fun calls
  20       iv(17)=maxfun
  21 * max num of iterations
  22       iv(18)=maxmin
  23 * controls output
  24       iv(19)=1
  25 * selects output unit
  26       iv(21)=out_minim
  27 * 1 means to print out result
  28       iv(22)=print_fin
  29 * 1 means to print out summary stats
  30       iv(23)=print_stat
  31 * 1 means to print initial x and d
  32       iv(24)=print_ini
  33 * min val for v(radfac) default is 0.1
  34       v(24)=0.1D0
  35 * max val for v(radfac) default is 4.0
  36       v(25)=2.0D0
  37 c     v(25)=4.0D0
  38 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
  39 * the sumsl default is 0.1
  40       v(26)=0.1D0
  41 * false conv if (act fnctn decrease) .lt. v(34)
  42 * the sumsl default is 100*machep
  43       v(34)=v(34)/100.0D0
  44 * absolute convergence
  45       if (tolf.eq.0.0D0) tolf=1.0D-6
  46       v(31)=tolf
  47 * relative convergence
  48       if (rtolf.eq.0.0D0) rtolf=1.0D-6
  49       v(32)=rtolf
  50 * controls initial step size
  51        v(35)=1.0D-2
  52 * large vals of d correspond to small components of step
  53       do i=1,n
  54         d(i)=1.0D0
  55       enddo
  56       return
  57       end
  58 c------------------------------------------------------------------------------
  59       subroutine minimize(n,x,f)
  60       implicit none
  61       include "DIMENSIONS"
  62       include "DIMENSIONS.ZSCOPT"
  63       integer liv,lv
  64       parameter (liv=60,lv=(77+max_paropt*(max_paropt+17)/2))
  65 *********************************************************************
  66 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
  67 * the calling subprogram.                                           *
  68 * when d(i)=1.0, then v(35) is the length of the initial step,      *
  69 * calculated in the usual pythagorean way.                          *
  70 * absolute convergence occurs when the function is within v(31) of  *
  71 * zero. unless you know the minimum value in advance, abs convg     *
  72 * is probably not useful.                                           *
  73 * relative convergence is when the model predicts that the function *
  74 * will decrease by less than v(32)*abs(fun).                        *
  75 *********************************************************************
  76 #ifdef MPI
  77       include "mpif.h"
  78       integer IERROR
  79       include "COMMON.MPI"
  80       integer status(MPI_STATUS_SIZE)
  81 #endif
  82       include "COMMON.WEIGHTS"
  83       include "COMMON.WEIGHTDER"
  84       include "COMMON.OPTIM"
  85       include "COMMON.CLASSES"
  86       include "COMMON.COMPAR"
  87       include "COMMON.MINPAR"
  88       include "COMMON.IOUNITS"
  89       include "COMMON.VMCPAR"
  90       include "COMMON.TIME1"
  91       include "COMMON.ENERGIES"
  92       include "COMMON.PROTNAME"
  93       include "COMMON.THERMAL"
  94       integer nn
  95       common /cc/ nn
  96       integer ilen
  97       external ilen
  98       integer iv(liv)
  99       double precision x(max_paropt),d(max_paropt),v(1:lv),g(max_paropt)
 100       real p(max_paropt+1,max_paropt),y(max_paropt+1),pb(max_paropt),
 101      &  pom(max_paropt),yb,y1,y2,pomi,ftol,temptr
 102       real funk
 103       external funk
 104       external targetfun,targetgrad,fdum
 105       integer uiparm,i,ii,j,k,nf,nfun,ngrad,iretcode,ianneal,maxanneal
 106       double precision urparm,penalty
 107       double precision f,f0,obf,obg(max_paropt),tcpu
 108       double precision Evarmax(maxprot),aux
 109       integer iscan,iprot,ibatch,niter
 110       integer n
 111       logical lprn
 112       lprn=.true.
 113       write (iout,*) "Minimizing with ",MINIMIZER
 114       IF (MINIMIZER.EQ."AMEBSA") THEN
 115         maxanneal=5
 116         nn=n
 117         do i=1,n
 118           do j=1,n+1
 119             p(j,i)=x(i)
 120           enddo
 121         enddo
 122         do j=1,n
 123           pom(j)=p(1,j)
 124         enddo
 125         y(1)=funk(pom)
 126         do i=2,n+1
 127           do j=1,n
 128             pom(j)=p(i,j)
 129           enddo
 130           pom(i-1)=p(i,i-1)-0.2
 131           pomi=pom(i-1)
 132           if (pom(i-1).lt.x_low(i-1)) pom(i-1)=x_low(i-1)
 133           y1=funk(pom)
 134           pom(i-1)=p(i,i-1)+0.2
 135           if (pom(i-1).gt.x_up(i-1)) pom(i-1)=x_up(i-1)
 136           y2=funk(pom)
 137           if (y2.lt.y1) then
 138             p(i,i-1)=pom(i-1)
 139             y(i)=y2
 140           else
 141             p(i,i-1)=pomi
 142             y(i)=y1
 143           endif
 144         enddo
 145         yb=1.0e20
 146         write (iout,*) "p and y"
 147         do i=1,n+1
 148           write (iout,'(20e15.5)') (p(i,j),j=1,n),y(i)
 149         enddo
 150         ftol=rtolf
 151         do ianneal=1,maxanneal
 152           iretcode=maxmin
 153           temptr=5.0/ianneal
 154           write (iout,*) "ianneal",ianneal," temptr",temptr
 155           call amebsa(p,y,max_paropt+1,max_paropt,n,pb,yb,ftol,funk,
 156      &      iretcode,temptr)
 157           write (iout,*) "Optimized P (PB)"
 158           do i=1,n
 159             write (iout,'(i5,f10.5)') i,pb(i)
 160           enddo
 161           write (iout,*) "Optimal function value",yb
 162           write(iout,*) "AMEBSA return code:",iretcode
 163         enddo
 164         do i=1,n
 165           x(i)=pb(i)
 166         enddo
 167         call targetfun(n,x,nf,f,uiparm,fdum)
 168         write (iout,*)
 169         write(iout,*) "Final function value:",f,f
 170         write(iout,*) "Final value of the object function:",obf
 171         write(iout,*) "Number of target function evaluations:",nfun
 172         write (iout,*)
 173         call x2w(n,x)
 174         do i=1, n
 175           xm(i)=x(i)
 176         enddo
 177         open (88,file=restartname,status="unknown")
 178         do i=1,n
 179           write (88,*) x(i)
 180         enddo
 181         close(88)
 182       ELSE IF (MINIMIZER.EQ."SUMSL") THEN
 183       nf=0
 184       nfun=0
 185       ngrad=0
 186       niter=0
 187       call minimize_init(n,iv,v,d)
 188 #ifdef CHECKGRAD
 189       write (iout,*) "Initial checkgrad:"
 190       call checkgrad(n,x)
 191 c      write (iout,*) "Second checkgrad"
 192 c      call checkgrad(n,x)
 193 #endif
 194       write (iout,*) "Calling SUMSL"
 195       iscan=0
 196       call targetfun(n,x,nf,f,uiparm,fdum)
 197       write (iout,*) "Initial function value",f
 198       f0=f
 199       f=-1.0d20
 200       do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 )
 201         f0=f
 202         call scan(n,x,f,.false.)
 203         iscan=iscan+1
 204       enddo
 205       PREVTIM=tcpu()
 206       f0=f
 207       f=-1.0d20
 208       if (maxscan.gt.0) then
 209         iscan=0
 210         do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 )
 211   101     continue
 212           cutoffviol=.false.
 213           call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm,
 214      &       urparm, fdum)
 215           write (iout,*) "Left SUMSL"
 216           penalty=v(10)
 217           iretcode=iv(1)
 218           nfun=iv(6)
 219           ngrad=iv(30)
 220           call targetfun(n,x,nf,f,uiparm,fdum)
 221           f0=f
 222           call scan(n,x,f,.false.)
 223           iscan=iscan+1
 224           if (iv(1).eq.11 .and. cutoffviol) then
 225             write (iout,*)
 226             write (iout,'(a)')
 227      &"Revising the list of conformations because of cut-off violation."
 228             do iprot=1,nprot
 229               write (iout,'(a,f7.1)')
 230      &          protname(iprot)(:ilen(protname(iprot))),enecut(iprot)
 231               do j=1,natlike(iprot)+2
 232                 write (iout,'(5(f7.1,f9.1,2x))')
 233      &           (1.0d0/(Rgas*betaT(k,ibatch,iprot)),
 234      &            Evar(j,k,iprot),k=1,nbeta(j,iprot))
 235               enddo
 236             enddo
 237 #ifdef MPI
 238             cutoffeval=.false.
 239             call func1(n,5,x)
 240 #else
 241             call make_list(.true.,.false.,n,x)
 242 #endif
 243             goto 101
 244           else if (iv(1).eq.11 .and. WhatsUp.eq.1) then
 245             write (iout,*) "Time up, terminating."
 246             goto 121
 247           endif
 248         enddo
 249       else
 250   102   continue
 251         cutoffviol=.false.
 252         cutoffeval=.false.
 253         write (iout,*) "MAXMIN",maxmin
 254         if (maxmin.gt.0)
 255      &  call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm,
 256      &   urparm, fdum)
 257         write (iout,*) "Left SUMSL"
 258         penalty=v(10)
 259         iretcode=iv(1)
 260         nfun=nfun+iv(6)
 261         ngrad=ngrad+iv(30)
 262         niter=niter+iv(31)
 263         write (iout,*) "iretcode",iv(1)," cutoffviol",cutoffviol
 264         if (iv(1).eq.11 .and. cutoffviol) then
 265           write (iout,'(a)')
 266      &"Revising the list of conformations because of cut-off violation."
 267           do iprot=1,nprot
 268             write (iout,'(a,f7.1)')
 269      &        protname(iprot)(:ilen(protname(iprot))),enecut(iprot)
 270             do j=1,natlike(iprot)+2
 271               write (iout,'(5(f7.1,f9.1,2x))')
 272      &         (1.0d0/(Rgas*betaT(k,ibatch,iprot)),
 273      &          Evar(j,k,iprot),k=1,nbeta(j,iprot))
 274             enddo
 275           enddo
 276           write (iout,'(a)')
 277           do iprot=1,nprot
 278             Evarmax(iprot)=Evar(1,1,iprot)*betaT(1,1,iprot)
 279             do ibatch=2,natlike(iprot)+2
 280               do k=1,nbeta(ibatch,iprot)
 281                 aux = Evar(k,ibatch,iprot)*betaT(k,ibatch,iprot)
 282                 if (aux.gt.Evarmax(iprot)) Evarmax(iprot)=aux
 283               enddo
 284             enddo
 285           enddo
 286           do iprot=1,nprot
 287             aux=dmin1(0.5d0*enecut_min(iprot)+1.25d0*Evarmax(iprot)
 288      &        ,enecut_max(iprot))
 289             enecut(iprot)=dmax1(aux,enecut_min(iprot))
 290           enddo
 291           write (iout,*) "Revised energy cutoffs"
 292           do iprot=1,nprot
 293             write (iout,*) "Protein",iprot,":",enecut(iprot)
 294           enddo
 295           call flush(iout)
 296 #ifdef MPI
 297 c
 298 c Send the order to the second master to revise the list of conformations.
 299 c
 300           cutoffeval=.false.
 301           call func1(n,5,x)
 302 #else
 303           call make_list(.true.,.false.,n,x)
 304 #endif
 305           call minimize_init(n,iv,v,d)
 306           iv(17)=iv(17)-nfun
 307           iv(18)=iv(18)-niter
 308           write (iout,*) "Diminishing IV(17) to",iv(17)
 309           write (iout,*) "Diminishing IV(18) to",iv(18)
 310           goto 102
 311         else if (iv(1).eq.11 .and. WhatsUp.eq.1) then
 312           write (iout,*) "Time up, terminating."
 313           goto 121
 314         endif
 315       endif
 316   121 continue
 317       call targetfun(n,x,nf,f,uiparm,fdum)
 318       write (iout,*)
 319       write(iout,*) "Final function value:",f,v(10)
 320       write(iout,*) "Final value of the object function:",obf
 321       write(iout,*) "Number of target function evaluations:",nfun
 322       write(iout,*) "Number of target gradient evaluations:",ngrad
 323       write(iout,*) "SUMSL return code:",iretcode
 324       write(*,*) "Final function value:",f,v(10)
 325       write(*,*) "Final value of the object function:",obf
 326       write(*,*) "Number of target function evaluations:",nfun
 327       write(*,*) "Number of target gradient evaluations:",ngrad
 328       write(*,*) "SUMSL return code:",iretcode
 329       write (iout,*)
 330       call x2w(n,x)
 331       do i=1, n
 332         xm(i)=x(i)
 333       enddo
 334       open (88,file=restartname,status="unknown")
 335       do i=1,n
 336         write (88,*) x(i)
 337       enddo
 338       close(88)
 339 #ifdef CHECKGRAD
 340       write (iout,*) "Final checkgrad:"
 341       call checkgrad(n,x)
 342 #endif
 343       ELSE
 344         write (iout,*) "Unknown minimizer."
 345       ENDIF
 346       return
 347       end
 348 c-----------------------------------------------------------------------------
 349       real function funk(p)
 350       implicit none
 351       include "DIMENSIONS"
 352       include "DIMENSIONS.ZSCOPT"
 353       real p(max_paropt)
 354       double precision x(max_paropt)
 355       double precision f
 356       integer uiparm
 357       double precision obf
 358       double precision ufparm
 359       include "COMMON.IOUNITS"
 360       integer n,nf
 361       common /cc/ n
 362       integer i
 363       external fdum
 364       nf=nf+1
 365       write (iout,*) "funk n",n," p",(p(i),i=1,n)
 366       do i=1,n
 367         x(i)=p(i)
 368       enddo
 369 c      write (iout,*) "funk n",n," x",(x(i),i=1,n)
 370       call targetfun(n,x,nf,f,uiparm,fdum)
 371       funk=f
 372       return
 373       end
 374 c-----------------------------------------------------------------------------
 375 c      logical function stopx(idum)
 376 c      integer idum
 377 c      stopx=.false.
 378 c      return
 379 c      end
 380 c-----------------------------------------------------------------------------
 381       subroutine checkgrad(n,x)
 382       implicit none
 383       include "DIMENSIONS"
 384       include "DIMENSIONS.ZSCOPT"
 385 #ifdef MPI
 386       include "mpif.h"
 387       include "COMMON.MPI"
 388       integer ierror
 389       integer status(MPI_STATUS_SIZE)
 390 #endif
 391       include "COMMON.IOUNITS"
 392       include "COMMON.VMCPAR"
 393       include "COMMON.OPTIM"
 394       include "COMMON.TIME1"
 395       include "COMMON.ENERGIES"
 396       include "COMMON.WEIGHTDER"
 397       include "COMMON.WEIGHTS"
 398       include "COMMON.CLASSES"
 399       include "COMMON.INTERACT"
 400       include "COMMON.NAMES"
 401       include "COMMON.VAR"
 402       include "COMMON.COMPAR"
 403       include "COMMON.TORSION"
 404       double precision energia(0:max_ene),epskl,eneps_num
 405       double precision f,fplus,fminus,fpot,fpotplus,xi,delta /5.0d-5/,
 406      & g(max_paropt),gnum(max_paropt),gfpot(max_paropt),gg(max_paropt),
 407      & rdum,zscder,pom,sumlikder1(maxvar,maxT,maxprot)
 408       integer n
 409       double precision x(n)
 410       integer uiparm,maxfun,i,ibatch,iprot,j,nf,k,kk,l,ii,iti,itj,ib,
 411      &  l1,l2,jj,inat
 412       double precision nuave_orig(maxdimnat,maxT,maxnatlike,maxprot)
 413 #ifdef EPSCHECK
 414       double precision enetb_orig1(maxstr_proc,max_ene,maxprot)
 415 #endif
 416       integer icant
 417       external icant
 418       double precision rdif
 419       double precision urparm,aux,fac,eoldsc
 420       double precision
 421      & rder_num,
 422      & rder_all(max_paropt,maxdimnat,maxT,maxnatlike,maxprot),
 423      & zvtot_orig(maxT,maxprot),sumlik_orig(maxT,maxprot)
 424       external fdum
 425       integer ind_shield /25/
 426
 427       write (iout,*) "Entered CHECKGRAD"
 428       call flush(iout)
 429
 430       call targetfun(n,x,nf,f,uiparm,fdum)
 431       call targetgrad(n,x,nf,g,uiparm,rdum,fdum)
 432       write (iout,*) "Initial function & gradient evaluation completed."
 433       call flush(iout)
 434 #ifdef EPSCHECK
 435       if (mod_side) then
 436       do iprot=1,nprot
 437         do i=1,ntot(iprot)
 438           jj = i2ii(i,iprot)
 439           if (jj.gt.0) then
 440           do j=1,n_ene
 441             enetb_orig1(jj,j,iprot)=enetb(jj,j,iprot)
 442           enddo
 443 c          write (iout,*) i,enetb_orig1(jj,1,iprot)
 444           endif
 445         enddo
 446       enddo
 447
 448       write (iout,*) "Derivatives in epsilons"
 449       DO iprot=1,nprot
 450
 451       call restore_molinfo(iprot)
 452       do i=1,ntot(iprot)
 453         jj = i2ii(i,iprot)
 454         if (jj.gt.0) then
 455         kk=0
 456 c         write (iout,*) "iprot",iprot," i",i
 457         call etotal(energia(0))
 458         eoldsc=energia(1)
 459         do k=1,ntyp
 460           do l=1,k
 461             epskl=eps(k,l)
 462             eps(k,l)=eps(k,l)+delta
 463             eps(l,k)=eps(k,l)
 464 c            write (iout,*) "***",k,l,epskl
 465             call prepare_sc
 466             kk=kk+1
 467             call restore_ccoords(iprot,i)
 468             call etotal(energia(0))
 469             eneps_num=(energia(1)-eoldsc)
 470      &        /delta
 471 c            write (iout,*) energia(1),enetb_orig1(jj,1,iprot)
 472             eps(k,l)=epskl
 473             eps(l,k)=epskl
 474 c            write (iout,*) "---",k,l,epskl
 475             write(iout,'(2a4,3e20.10)') restyp(k),restyp(l),
 476      &        eneps(kk,i,iprot),eneps_num,
 477      &        (eneps_num-eneps(kk,i,iprot))*100
 478           enddo
 479         enddo
 480         endif
 481       enddo
 482
 483       ENDDO
 484
 485       call prepare_sc
 486
 487       endif
 488 #endif
 489       write (iout,*) "calling func1 0"
 490       call flush(iout)
 491       call func1(n,1,x)
 492 #ifdef NEWCORR
 493 c AL 2/10/2018 Add PMF gradient
 494       if (mod_fourier(nloctyp).gt.0) then
 495         fpot = rdif(n,x,gfpot,2)
 496       endif
 497 #endif
 498       do iprot=1,nprot
 499         do inat=1,natlike(iprot)
 500           do ib=1,nbeta(inat+2,iprot)
 501             do i=1,natdim(inat,iprot)
 502               nuave_orig(i,ib,inat,iprot)=nuave(i,ib,inat,iprot)
 503             enddo
 504           enddo
 505         enddo
 506       enddo
 507       do iprot=1,nprot
 508         do ib=1,nbeta(1,iprot)
 509           sumlik_orig(ib,iprot)=sumlik(ib,iprot)
 510         enddo
 511       enddo
 512       do iprot=1,nprot
 513         do ib=1,nbeta(2,iprot)
 514           zvtot_orig(ib,iprot)=zvtot(ib,iprot)
 515         enddo
 516       enddo
 517       write(iout,*) "exit func1 0"
 518       do iprot=1,nprot
 519         do ib=1,nbeta(2,iprot)
 520           do i=1,n
 521             zvdertot(i,ib,iprot)=0.0d0
 522           enddo
 523         enddo
 524       enddo
 525       call zderiv
 526       ii=0
 527       do k=1,n_ene
 528         if (imask(k).gt.0 .and. k.ne.ind_shield) ii=ii+1
 529       enddo
 530       do iprot=1,nprot
 531         do ib=1,nbeta(1,iprot)
 532           do k=1,ii
 533             sumlikder1(k,ib,iprot)=sumlikder(k,ib,iprot)
 534 c            write (iout,*) "ib",ib," k",k," sumlikder",
 535 c     &            sumlikder1(k,ib,iprot)
 536           enddo
 537         enddo
 538       enddo
 539 #ifdef NEWCORR
 540 c AL 2/10/2018 Add PMF gradient
 541       if (mod_fourier(nloctyp).gt.0 .and. torsion_pmf) then
 542         do i=1,nloctyp*(2*nloctyp-1)
 543           ii=ii+1
 544           sumlikder1(ii,ib,iprot)=0.0d0
 545         enddo
 546       endif
 547 #endif
 548       do iprot=1,nprot
 549 c Maximum likelihood terms
 550         do ib=1,nbeta(1,iprot)
 551           kk=ii
 552           do i=1,nsingle_sc
 553             kk=kk+1
 554             iti=ityp_ssc(i)
 555             sumlikder1(kk,ib,iprot)=
 556      &       sumlikeps(icant(iti,iti),ib,iprot)
 557             do j=1,ntyp
 558               sumlikder1(kk,ib,iprot)=
 559      &         sumlikder1(kk,ib,iprot)+
 560      &         sumlikeps(icant(iti,itj),ib,iprot)
 561             enddo
 562           enddo
 563           do i=1,npair_sc
 564             kk=kk+1
 565             iti=ityp_psc(1,i)
 566             itj=ityp_psc(2,i)
 567             sumlikder1(kk,ib,iprot)=
 568      &        sumlikeps(icant(iti,itj),ib,iprot)
 569 c            write (iout,*) "kk",kk," iprot",iprot," ib",ib,
 570 c     &       " iti",iti," itj",itj,
 571 c     &       " ind",icant(iti,itj)," sumlikeps",
 572 c     &       sumlikeps(icant(iti,itj),ib,iprot)
 573           enddo
 574           do i=1,ntor_var
 575             kk=kk+1
 576             sumlikder1(kk,ib,iprot)=sumliktor(i,ib,iprot)
 577           enddo
 578           do i=1,nang_var
 579             kk=kk+1
 580             sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot)
 581           enddo
 582         enddo
 583 c Heat capacity terms
 584         do ib=1,nbeta(2,iprot)
 585           kk=ii
 586           do i=1,nsingle_sc
 587             kk=kk+1
 588             iti=ityp_ssc(i)
 589             zvdertot(kk,ib,iprot)=
 590      &       zvepstot(icant(iti,iti),ib,iprot)
 591             do j=1,ntyp
 592               zvdertot(kk,ib,iprot)=
 593      &         zvdertot(kk,ib,iprot)+
 594      &         zvepstot(icant(iti,j),ib,iprot)
 595             enddo
 596           enddo
 597           do i=1,npair_sc
 598             kk=kk+1
 599             iti=ityp_psc(1,i)
 600             itj=ityp_psc(2,i)
 601             zvdertot(kk,ib,iprot)=
 602      &        zvepstot(icant(iti,itj),ib,iprot)
 603           enddo
 604           do i=1,ntor_var
 605             kk=kk+1
 606             zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
 607           enddo
 608           do i=1,nang_var
 609             kk=kk+1
 610             zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot)
 611           enddo
 612         enddo
 613       enddo
 614 c Molecular properties
 615       do iprot=1,nprot
 616         do inat=1,natlike(iprot)
 617           do ib=1,nbeta(inat+2,iprot)
 618             call propderiv(ib,inat,iprot)
 619             do k=1,natdim(inat,iprot)
 620               do kk=1,ii
 621                 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
 622               enddo
 623               kk=ii
 624               do i=1,nsingle_sc
 625                 kk=kk+1
 626                 iti=ityp_ssc(i)
 627                 rder_all(kk,k,ib,inat,iprot)=
 628      &           reps(icant(iti,iti),k)
 629                 do j=1,ntyp
 630                   rder_all(kk,k,ib,inat,iprot)=
 631      &             rder_all(kk,k,ib,inat,iprot)+
 632      &             reps(icant(iti,j),k)
 633                 enddo
 634               enddo
 635               do i=1,npair_sc
 636                 kk=kk+1
 637                 iti=ityp_psc(1,i)
 638                 itj=ityp_psc(2,i)
 639                 rder_all(kk,k,ib,inat,iprot)=
 640      &           reps(icant(iti,itj),k)
 641               enddo
 642               kk=ii
 643               do i=1,ntor_var
 644                 kk=kk+1
 645                 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
 646               enddo
 647               do i=1,nang_var
 648                 kk=kk+1
 649                 rder_all(kk,k,ib,inat,iprot)=rang(i,k)
 650               enddo
 651             enddo
 652           enddo
 653         enddo
 654       enddo
 655 c Calculate numerical derivatives and compare with analytical derivatives
 656       do i=1,kk
 657         xi=x(i)
 658         x(i)=xi+xi*delta
 659 c        write (iout,*) (x(k),k=1,kk)
 660         call func1(n,1,x)
 661 c        write (iout,*) "after func1"
 662 c        call flush(iout)
 663 #ifdef NEWCORR
 664 c AL 2/10/2018 Add PMF gradient
 665         if (mod_fourier(nloctyp).gt.0) then
 666           fpotplus = rdif(n,x,gfpot,1)
 667           write (iout,*) "after rdif"
 668           call flush(iout)
 669         endif
 670 #endif
 671         write (iout,*) "Variable",i
 672 c Maximum likelihood
 673         do iprot=1,nprot
 674           write (iout,*) "Derivatives of maximum likelihood"
 675           do ib=1,nbeta(1,iprot)
 676 c            write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
 677             pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta)
 678             write (iout,'(2i5,3e14.5)') iprot,ib,
 679      &        sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
 680      &        /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
 681 c            write (iout,*) sumlik(ib,iprot),
 682 c     &        sumlik_orig(ib,iprot)
 683           enddo
 684         enddo
 685 c Heat capacity
 686         do iprot=1,nprot
 687           write (iout,*) "Derivatives of heat capacity"
 688           do ib=1,nbeta(2,iprot)
 689             pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta)
 690             write (iout,'(2i5,3e14.5)') iprot,ib,
 691      &        zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
 692      &        /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
 693 c            write (iout,*) zvtot(ib,ibatch,iprot),
 694 c     &        zvtot_orig(ib,ibatch,iprot)
 695           enddo
 696         enddo
 697 c Numerical derivatives of natlike properties
 698         write (iout,*) "Derivatives of molecular properties"
 699         do iprot=1,nprot
 700           do inat=1,natlike(iprot)
 701             do ib=1,nbeta(inat+2,iprot)
 702               do k=1,natdim(inat,iprot)
 703                 rder_num=(nuave(k,ib,inat,iprot)
 704      &                   -nuave_orig(k,ib,inat,iprot))/(xi*delta)
 705                 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
 706      &             rder_all(kk,k,ib,inat,iprot),
 707      &             rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
 708      &             (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
 709               enddo
 710             enddo
 711           enddo
 712         enddo
 713 #ifdef NEWCORR
 714 c AL 2/10/2018 Add PMF gradient
 715         if (mod_fourier(nloctyp).gt.0) then
 716           write(iout,*)"Derivative of fpot",(-fpotplus+fpot)/(xi*delta),
 717      &      gfpot(i),(gfpot(i)-(-fpotplus+fpot)/(xi*delta))/gfpot(i)
 718         endif
 719 #endif
 720         x(i)=xi
 721       enddo
 722       write (iout,*)
 723       do i=1,n
 724         xi=x(i)
 725         x(i)=xi+xi*delta
 726 c        write (iout,*) "i",i," xplus ",x(i),xi
 727         call targetfun(n,x,nf,fplus,uiparm,fdum)
 728         x(i)=xi-xi*delta
 729 c        write (iout,*) "i",i," xminus",x(i),xi
 730         call targetfun(n,x,nf,fminus,uiparm,fdum)
 731         x(i)=xi
 732 c        write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus,
 733 c     &   "obfplus",obfplus," obfminus",obfminus
 734         gnum(i)=0.5d0*(fplus-fminus)/(xi*delta)
 735       enddo
 736       write (iout,*) "Comparison of analytical and numerical gradient"
 737       do i=1,n
 738         write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
 739      &    1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
 740       enddo
 741       write (iout,*)
 742       return
 743       end
 744 c-----------------------------------------------------------------------
 745       subroutine prepare_sc
 746       implicit none
 747       include "DIMENSIONS"
 748       include "DIMENSIONS.ZSCOPT"
 749       include "COMMON.NAMES"
 750       include "COMMON.WEIGHTS"
 751       include "COMMON.INTERACT"
 752       include "COMMON.ENERGIES"
 753       include "COMMON.FFIELD"
 754       include "COMMON.TORSION"
 755       include "COMMON.IOUNITS"
 756       logical lprint
 757       integer i,j,ii,nvarr,iti,itj
 758       double precision rri
 759       double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
 760      &  ratsig1,ratsig2,rsum_max,dwa16,r_augm
 761
 762       lprint=.false.
 763
 764         if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
 765      &    potname(ipot),', exponents are ',expon,2*expon
 766 C-----------------------------------------------------------------------
 767 C Calculate the "working" parameters of SC interactions.
 768         dwa16=2.0d0**(1.0d0/6.0d0)
 769         do i=1,ntyp
 770           do j=i,ntyp
 771             sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 772             sigma(j,i)=sigma(i,j)
 773             rs0(i,j)=dwa16*sigma(i,j)
 774             rs0(j,i)=rs0(i,j)
 775           enddo
 776         enddo
 777         if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 778      &   'Working parameters of the SC interactions:',
 779      &   '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 780      &   '  chi1   ','   chi2   '
 781         do i=1,ntyp
 782           do j=i,ntyp
 783             epsij=eps(i,j)
 784             if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 785               rrij=sigma(i,j)
 786             else
 787               rrij=rr0(i)+rr0(j)
 788             endif
 789             r0(i,j)=rrij
 790             r0(j,i)=rrij
 791             rrij=rrij**expon
 792             epsij=eps(i,j)
 793             sigeps=dsign(1.0D0,epsij)
 794             epsij=dabs(epsij)
 795             aa(i,j)=epsij*rrij*rrij
 796             bb(i,j)=-sigeps*epsij*rrij
 797             aa(j,i)=aa(i,j)
 798             bb(j,i)=bb(i,j)
 799             if (ipot.gt.2) then
 800               sigt1sq=sigma0(i)**2
 801               sigt2sq=sigma0(j)**2
 802               sigii1=sigii(i)
 803               sigii2=sigii(j)
 804               ratsig1=sigt2sq/sigt1sq
 805               ratsig2=1.0D0/ratsig1
 806               chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 807               if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 808               rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 809             else
 810               rsum_max=sigma(i,j)
 811             endif
 812             sigmaii(i,j)=rsum_max
 813             sigmaii(j,i)=rsum_max
 814             if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
 815      &      then
 816               r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 817               augm(i,j)=epsij*r_augm**(2*expon)
 818 c             augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 819               augm(j,i)=augm(i,j)
 820             else
 821               augm(i,j)=0.0D0
 822               augm(j,i)=0.0D0
 823             endif
 824             if (lprint) then
 825               write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 826      &        restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 827      &        sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 828             endif
 829           enddo
 830         enddo
 831       return
 832       end