source/maxlik/src-Fmatch/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-7/,
 406      & g(max_paropt),gnum(max_paropt),gfpot(max_paropt),gg(max_paropt),
 407      & gforce(max_paropt),rdum,zscder,pom,fchi,fchiplus,
 408      & sumlikder1(maxvar,maxT,maxprot)
 409       double precision chisquare_force
 410       integer n
 411       double precision x(n)
 412       integer uiparm,maxfun,i,ibatch,iprot,j,nf,k,kk,l,ii,iti,itj,ib,
 413      &  l1,l2,jj,inat
 414       double precision nuave_orig(maxdimnat,maxT,maxnatlike,maxprot)
 415 #ifdef EPSCHECK
 416       double precision enetb_orig1(maxstr_proc,max_ene,maxprot)
 417 #endif
 418       integer icant
 419       external icant
 420       double precision rdif
 421       double precision urparm,aux,fac,eoldsc
 422       double precision
 423      & rder_num,
 424      & rder_all(max_paropt,maxdimnat,maxT,maxnatlike,maxprot),
 425      & zvtot_orig(maxT,maxprot),sumlik_orig(maxT,maxprot)
 426       external fdum
 427       integer ind_shield /25/
 428
 429 c      write (iout,*) "Entered CHECKGRAD"
 430       call flush(iout)
 431
 432       call targetfun(n,x,nf,f,uiparm,fdum)
 433 c      write (iout,*) "f",f
 434       call targetgrad(n,x,nf,g,uiparm,rdum,fdum)
 435 c      write (iout,*) "g",g(:n)
 436       write (iout,*) "Initial function & gradient evaluation completed."
 437       call flush(iout)
 438 #ifdef EPSCHECK
 439       if (mod_side) then
 440       do iprot=1,nprot
 441         do i=1,ntot(iprot)
 442           jj = i2ii(i,iprot)
 443           if (jj.gt.0) then
 444           do j=1,n_ene
 445             enetb_orig1(jj,j,iprot)=enetb(jj,j,iprot)
 446           enddo
 447 c          write (iout,*) i,enetb_orig1(jj,1,iprot)
 448           endif
 449         enddo
 450       enddo
 451
 452       write (iout,*) "Derivatives in epsilons"
 453       DO iprot=1,nprot
 454
 455       call restore_molinfo(iprot)
 456       do i=1,ntot(iprot)
 457         jj = i2ii(i,iprot)
 458         if (jj.gt.0) then
 459         kk=0
 460 c         write (iout,*) "iprot",iprot," i",i
 461         call etotal(energia(0))
 462         eoldsc=energia(1)
 463         do k=1,ntyp
 464           do l=1,k
 465             epskl=eps(k,l)
 466             eps(k,l)=eps(k,l)+delta
 467             eps(l,k)=eps(k,l)
 468 c            write (iout,*) "***",k,l,epskl
 469             call prepare_sc
 470             kk=kk+1
 471             call restore_ccoords(iprot,i)
 472             call etotal(energia(0))
 473             eneps_num=(energia(1)-eoldsc)
 474      &        /delta
 475 c            write (iout,*) energia(1),enetb_orig1(jj,1,iprot)
 476             eps(k,l)=epskl
 477             eps(l,k)=epskl
 478 c            write (iout,*) "---",k,l,epskl
 479             write(iout,'(2a4,3e20.10)') restyp(k),restyp(l),
 480      &        eneps(kk,i,iprot),eneps_num,
 481      &        (eneps_num-eneps(kk,i,iprot))*100
 482           enddo
 483         enddo
 484         endif
 485       enddo
 486
 487       ENDDO
 488
 489       call prepare_sc
 490
 491       endif
 492 #endif
 493       write (iout,*) "calling func1 0"
 494       call flush(iout)
 495       call func1(n,1,x)
 496 #ifdef NEWCORR
 497 c AL 2/10/2018 Add PMF gradient
 498       if (mod_fourier(nloctyp).gt.0) then
 499         fpot = rdif(n,x,gfpot,2)
 500       endif
 501 #endif
 502 c      write (iout,*) "-----------------------------------"
 503 c      write (iout,*) "Computing force chisquare"
 504       fchi = chisquare_force(n,x,gforce,2)
 505 c      write (iout,*) "After first call to chisquare_force"
 506 c      write (iout,*) "fchi",fchi
 507 c      write (iout,*) "gforce",gforce(:n)
 508 c      write (iout,*) "-----------------------------------"
 509       do iprot=1,nprot
 510         do inat=1,natlike(iprot)
 511           do ib=1,nbeta(inat+2,iprot)
 512             do i=1,natdim(inat,iprot)
 513               nuave_orig(i,ib,inat,iprot)=nuave(i,ib,inat,iprot)
 514             enddo
 515           enddo
 516         enddo
 517       enddo
 518       do iprot=1,nprot
 519         do ib=1,nbeta(1,iprot)
 520           sumlik_orig(ib,iprot)=sumlik(ib,iprot)
 521         enddo
 522       enddo
 523       do iprot=1,nprot
 524         do ib=1,nbeta(2,iprot)
 525           zvtot_orig(ib,iprot)=zvtot(ib,iprot)
 526         enddo
 527       enddo
 528 c      write(iout,*) "exit func1 0"
 529       do iprot=1,nprot
 530         do ib=1,nbeta(2,iprot)
 531           do i=1,n
 532             zvdertot(i,ib,iprot)=0.0d0
 533           enddo
 534         enddo
 535       enddo
 536       call zderiv
 537       ii=0
 538       do k=1,n_ene
 539         if (imask(k).gt.0 .and. k.ne.ind_shield) ii=ii+1
 540       enddo
 541       do iprot=1,nprot
 542         do ib=1,nbeta(1,iprot)
 543           do k=1,ii
 544             sumlikder1(k,ib,iprot)=sumlikder(k,ib,iprot)
 545 c            write (iout,*) "ib",ib," k",k," sumlikder",
 546 c     &            sumlikder1(k,ib,iprot)
 547           enddo
 548         enddo
 549       enddo
 550 #ifdef NEWCORR
 551 c AL 2/10/2018 Add PMF gradient
 552       if (mod_fourier(nloctyp).gt.0 .and. torsion_pmf) then
 553         do i=1,nloctyp*(2*nloctyp-1)
 554           ii=ii+1
 555           sumlikder1(ii,ib,iprot)=0.0d0
 556         enddo
 557       endif
 558 #endif
 559       kk=ii
 560       do iprot=1,nprot
 561 c Maximum likelihood terms
 562         do ib=1,nbeta(1,iprot)
 563           kk=ii
 564           do i=1,nsingle_sc
 565             kk=kk+1
 566             iti=ityp_ssc(i)
 567             sumlikder1(kk,ib,iprot)=
 568      &       sumlikeps(icant(iti,iti),ib,iprot)
 569             do j=1,ntyp
 570               sumlikder1(kk,ib,iprot)=
 571      &         sumlikder1(kk,ib,iprot)+
 572      &         sumlikeps(icant(iti,itj),ib,iprot)
 573             enddo
 574           enddo
 575           do i=1,npair_sc
 576             kk=kk+1
 577             iti=ityp_psc(1,i)
 578             itj=ityp_psc(2,i)
 579             sumlikder1(kk,ib,iprot)=
 580      &        sumlikeps(icant(iti,itj),ib,iprot)
 581 c            write (iout,*) "kk",kk," iprot",iprot," ib",ib,
 582 c     &       " iti",iti," itj",itj,
 583 c     &       " ind",icant(iti,itj)," sumlikeps",
 584 c     &       sumlikeps(icant(iti,itj),ib,iprot)
 585           enddo
 586           do i=1,ntor_var
 587             kk=kk+1
 588             sumlikder1(kk,ib,iprot)=sumliktor(i,ib,iprot)
 589           enddo
 590           do i=1,nang_var
 591             kk=kk+1
 592             sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot)
 593           enddo
 594         enddo
 595 c Heat capacity terms
 596         do ib=1,nbeta(2,iprot)
 597           kk=ii
 598           do i=1,nsingle_sc
 599             kk=kk+1
 600             iti=ityp_ssc(i)
 601             zvdertot(kk,ib,iprot)=
 602      &       zvepstot(icant(iti,iti),ib,iprot)
 603             do j=1,ntyp
 604               zvdertot(kk,ib,iprot)=
 605      &         zvdertot(kk,ib,iprot)+
 606      &         zvepstot(icant(iti,j),ib,iprot)
 607             enddo
 608           enddo
 609           do i=1,npair_sc
 610             kk=kk+1
 611             iti=ityp_psc(1,i)
 612             itj=ityp_psc(2,i)
 613             zvdertot(kk,ib,iprot)=
 614      &        zvepstot(icant(iti,itj),ib,iprot)
 615           enddo
 616           do i=1,ntor_var
 617             kk=kk+1
 618             zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
 619           enddo
 620           do i=1,nang_var
 621             kk=kk+1
 622             zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot)
 623           enddo
 624         enddo
 625       enddo
 626 c Molecular properties
 627       do iprot=1,nprot
 628         do inat=1,natlike(iprot)
 629           do ib=1,nbeta(inat+2,iprot)
 630             call propderiv(ib,inat,iprot)
 631             do k=1,natdim(inat,iprot)
 632               do kk=1,ii
 633                 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
 634               enddo
 635               kk=ii
 636               do i=1,nsingle_sc
 637                 kk=kk+1
 638                 iti=ityp_ssc(i)
 639                 rder_all(kk,k,ib,inat,iprot)=
 640      &           reps(icant(iti,iti),k)
 641                 do j=1,ntyp
 642                   rder_all(kk,k,ib,inat,iprot)=
 643      &             rder_all(kk,k,ib,inat,iprot)+
 644      &             reps(icant(iti,j),k)
 645                 enddo
 646               enddo
 647               do i=1,npair_sc
 648                 kk=kk+1
 649                 iti=ityp_psc(1,i)
 650                 itj=ityp_psc(2,i)
 651                 rder_all(kk,k,ib,inat,iprot)=
 652      &           reps(icant(iti,itj),k)
 653               enddo
 654               kk=ii
 655               do i=1,ntor_var
 656                 kk=kk+1
 657                 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
 658               enddo
 659               do i=1,nang_var
 660                 kk=kk+1
 661                 rder_all(kk,k,ib,inat,iprot)=rang(i,k)
 662               enddo
 663             enddo
 664           enddo
 665         enddo
 666       enddo
 667 c Calculate numerical derivatives and compare with analytical derivatives
 668       do i=1,kk
 669 c        write (iout,*) "gforce1:",gforce(:n)
 670         xi=x(i)
 671         x(i)=xi+xi*delta
 672 c        write (iout,*) (x(k),k=1,kk)
 673         call func1(n,1,x)
 674
 675 c        write (iout,*) "after func1"
 676 c        call flush(iout)
 677 #ifdef NEWCORR
 678 c AL 2/10/2018 Add PMF gradient
 679         if (mod_fourier(nloctyp).gt.0) then
 680           fpotplus = rdif(n,x,gfpot,1)
 681 c          write (iout,*) "after rdif"
 682 c          call flush(iout)
 683         endif
 684 #endif
 685         fchiplus = chisquare_force(n,x,gg,1)
 686 c        write (iout,*) "gforce2:",gforce(:n)
 687 c        write (iout,*) "Variable",i," fchiplus",fchiplus
 688 c Maximum likelihood
 689         do iprot=1,nprot
 690           write (iout,*) "Derivatives of maximum likelihood"
 691           do ib=1,nbeta(1,iprot)
 692 c            write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
 693             pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta)
 694             write (iout,'(2i5,3e14.5)') iprot,ib,
 695      &        sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
 696      &        /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
 697 c            write (iout,*) sumlik(ib,iprot),
 698 c     &        sumlik_orig(ib,iprot)
 699           enddo
 700         enddo
 701 c Heat capacity
 702         do iprot=1,nprot
 703           write (iout,*) "Derivatives of heat capacity"
 704           do ib=1,nbeta(2,iprot)
 705             pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta)
 706             write (iout,'(2i5,3e14.5)') iprot,ib,
 707      &        zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
 708      &        /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
 709 c            write (iout,*) zvtot(ib,ibatch,iprot),
 710 c     &        zvtot_orig(ib,ibatch,iprot)
 711           enddo
 712         enddo
 713 c Numerical derivatives of natlike properties
 714         write (iout,*) "Derivatives of molecular properties"
 715         do iprot=1,nprot
 716           do inat=1,natlike(iprot)
 717             do ib=1,nbeta(inat+2,iprot)
 718               do k=1,natdim(inat,iprot)
 719                 rder_num=(nuave(k,ib,inat,iprot)
 720      &                   -nuave_orig(k,ib,inat,iprot))/(xi*delta)
 721                 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
 722      &             rder_all(kk,k,ib,inat,iprot),
 723      &             rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
 724      &             (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
 725               enddo
 726             enddo
 727           enddo
 728         enddo
 729 #ifdef NEWCORR
 730 c AL 2/10/2018 Add PMF gradient
 731         if (mod_fourier(nloctyp).gt.0) then
 732           write(iout,*)"Derivative of fpot",(-fpotplus+fpot)/(xi*delta),
 733      &      gfpot(i),(gfpot(i)-(-fpotplus+fpot)/(xi*delta))/gfpot(i)
 734         endif
 735 #endif
 736 c        write (iout,*) "fchiplus",fchiplus," fchi",fchi," xi",xi,
 737 c     &    " delta",delta
 738         write (iout,*)"Derivatives of fchi",(-fchiplus+fchi)/(xi*delta),
 739      &      gforce(i),(gforce(i)-(-fchiplus+fchi)/(xi*delta))/gforce(i)
 740         x(i)=xi
 741       enddo
 742       write (iout,*)
 743       do i=1,n
 744         xi=x(i)
 745         x(i)=xi+xi*delta
 746 c        write (iout,*) "i",i," xplus ",x(i),xi
 747         call targetfun(n,x,nf,fplus,uiparm,fdum)
 748         x(i)=xi-xi*delta
 749 c        write (iout,*) "i",i," xminus",x(i),xi
 750         call targetfun(n,x,nf,fminus,uiparm,fdum)
 751         x(i)=xi
 752         write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus
 753         gnum(i)=0.5d0*(fplus-fminus)/(xi*delta)
 754       enddo
 755       write (iout,*) "Comparison of analytical and numerical gradient"
 756       do i=1,n
 757         write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
 758      &    1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
 759       enddo
 760       write (iout,*)
 761       return
 762       end
 763 c-----------------------------------------------------------------------
 764       subroutine prepare_sc
 765       implicit none
 766       include "DIMENSIONS"
 767       include "DIMENSIONS.ZSCOPT"
 768       include "COMMON.NAMES"
 769       include "COMMON.WEIGHTS"
 770       include "COMMON.INTERACT"
 771       include "COMMON.ENERGIES"
 772       include "COMMON.FFIELD"
 773       include "COMMON.TORSION"
 774       include "COMMON.IOUNITS"
 775       logical lprint
 776       integer i,j,ii,nvarr,iti,itj
 777       double precision rri
 778       double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
 779      &  ratsig1,ratsig2,rsum_max,dwa16,r_augm
 780
 781       lprint=.false.
 782
 783         if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
 784      &    potname(ipot),', exponents are ',expon,2*expon
 785 C-----------------------------------------------------------------------
 786 C Calculate the "working" parameters of SC interactions.
 787         dwa16=2.0d0**(1.0d0/6.0d0)
 788         do i=1,ntyp
 789           do j=i,ntyp
 790             sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 791             sigma(j,i)=sigma(i,j)
 792             rs0(i,j)=dwa16*sigma(i,j)
 793             rs0(j,i)=rs0(i,j)
 794           enddo
 795         enddo
 796         if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 797      &   'Working parameters of the SC interactions:',
 798      &   '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 799      &   '  chi1   ','   chi2   '
 800         do i=1,ntyp
 801           do j=i,ntyp
 802             epsij=eps(i,j)
 803             if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 804               rrij=sigma(i,j)
 805             else
 806               rrij=rr0(i)+rr0(j)
 807             endif
 808             r0(i,j)=rrij
 809             r0(j,i)=rrij
 810             rrij=rrij**expon
 811             epsij=eps(i,j)
 812             sigeps=dsign(1.0D0,epsij)
 813             epsij=dabs(epsij)
 814             aa(i,j)=epsij*rrij*rrij
 815             bb(i,j)=-sigeps*epsij*rrij
 816             aa(j,i)=aa(i,j)
 817             bb(j,i)=bb(i,j)
 818             if (ipot.gt.2) then
 819               sigt1sq=sigma0(i)**2
 820               sigt2sq=sigma0(j)**2
 821               sigii1=sigii(i)
 822               sigii2=sigii(j)
 823               ratsig1=sigt2sq/sigt1sq
 824               ratsig2=1.0D0/ratsig1
 825               chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 826               if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 827               rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 828             else
 829               rsum_max=sigma(i,j)
 830             endif
 831             sigmaii(i,j)=rsum_max
 832             sigmaii(j,i)=rsum_max
 833             if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
 834      &      then
 835               r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 836               augm(i,j)=epsij*r_augm**(2*expon)
 837 c             augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 838               augm(j,i)=augm(i,j)
 839             else
 840               augm(i,j)=0.0D0
 841               augm(j,i)=0.0D0
 842             endif
 843             if (lprint) then
 844               write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 845      &        restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 846      &        sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 847             endif
 848           enddo
 849         enddo
 850       return
 851       end