source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/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,fpdbplus,xi,
 406      & delta /5.0d-7/,
 407      & g(max_paropt),gnum(max_paropt),gfpot(max_paropt),
 408      & gfpdb(max_paropt),gg(max_paropt),
 409      & rdum,zscder,pom,sumlikder1(maxvar,maxT,maxprot)
 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,maxlik_pdb
 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       write (iout,*) "Entered CHECKGRAD"
 430       call flush(iout)
 431
 432       call targetfun(n,x,nf,f,uiparm,fdum)
 433       write (iout,*) "Exitted targetfun"
 434       call flush(iout)
 435       call targetgrad(n,x,nf,g,uiparm,rdum,fdum)
 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 c      if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 499       if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 500      &  fpot = rdif(n,x,gfpot,2)
 501       if (pdbpmf) fpdb = maxlik_pdb(n,x,gfpdb,2)
 502 c      endif
 503 #endif
 504       do iprot=1,nprot
 505         do inat=1,natlike(iprot)
 506           do ib=1,nbeta(inat+2,iprot)
 507             do i=1,natdim(inat,iprot)
 508               nuave_orig(i,ib,inat,iprot)=nuave(i,ib,inat,iprot)
 509             enddo
 510           enddo
 511         enddo
 512       enddo
 513       do iprot=1,nprot
 514         do ib=1,nbeta(1,iprot)
 515           sumlik_orig(ib,iprot)=sumlik(ib,iprot)
 516         enddo
 517       enddo
 518       do iprot=1,nprot
 519         do ib=1,nbeta(2,iprot)
 520           zvtot_orig(ib,iprot)=zvtot(ib,iprot)
 521         enddo
 522       enddo
 523       write(iout,*) "exit func1 0"
 524       do iprot=1,nprot
 525         do ib=1,nbeta(2,iprot)
 526           do i=1,n
 527             zvdertot(i,ib,iprot)=0.0d0
 528           enddo
 529         enddo
 530       enddo
 531       call zderiv
 532       ii=0
 533       do k=1,n_ene
 534         if (imask(k).gt.0 .and. k.ne.ind_shield) ii=ii+1
 535       enddo
 536       do iprot=1,nprot
 537         do ib=1,nbeta(1,iprot)
 538           do k=1,ii
 539             sumlikder1(k,ib,iprot)=sumlikder(k,ib,iprot)
 540 c            write (iout,*) "ib",ib," k",k," sumlikder",
 541 c     &            sumlikder1(k,ib,iprot)
 542           enddo
 543         enddo
 544       enddo
 545 #ifdef NEWCORR
 546 c AL 2/10/2018 Add PMF gradient
 547       if ((mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)
 548      &    .and.(torsion_pmf.or.turn_pmf.or.eello_pmf)) then
 549         do i=1,nloctyp*(2*nloctyp-1)
 550           ii=ii+1
 551           sumlikder1(ii,ib,iprot)=0.0d0
 552         enddo
 553       endif
 554 #endif
 555       do iprot=1,nprot
 556 c Maximum likelihood terms
 557         do ib=1,nbeta(1,iprot)
 558           kk=ii
 559           do i=1,nsingle_sc
 560             kk=kk+1
 561             iti=ityp_ssc(i)
 562             sumlikder1(kk,ib,iprot)=
 563      &       sumlikeps(icant(iti,iti),ib,iprot)
 564             do j=1,ntyp
 565               sumlikder1(kk,ib,iprot)=
 566      &         sumlikder1(kk,ib,iprot)+
 567      &         sumlikeps(icant(iti,itj),ib,iprot)
 568             enddo
 569           enddo
 570           do i=1,npair_sc
 571             kk=kk+1
 572             iti=ityp_psc(1,i)
 573             itj=ityp_psc(2,i)
 574             sumlikder1(kk,ib,iprot)=
 575      &        sumlikeps(icant(iti,itj),ib,iprot)
 576 c            write (iout,*) "kk",kk," iprot",iprot," ib",ib,
 577 c     &       " iti",iti," itj",itj,
 578 c     &       " ind",icant(iti,itj)," sumlikeps",
 579 c     &       sumlikeps(icant(iti,itj),ib,iprot)
 580           enddo
 581           do i=1,ntor_var
 582             kk=kk+1
 583             sumlikder1(kk,ib,iprot)=sumliktor(i,ib,iprot)
 584           enddo
 585           do i=1,nang_var
 586             kk=kk+1
 587             sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot)
 588           enddo
 589         enddo
 590 c Heat capacity terms
 591         do ib=1,nbeta(2,iprot)
 592           kk=ii
 593           do i=1,nsingle_sc
 594             kk=kk+1
 595             iti=ityp_ssc(i)
 596             zvdertot(kk,ib,iprot)=
 597      &       zvepstot(icant(iti,iti),ib,iprot)
 598             do j=1,ntyp
 599               zvdertot(kk,ib,iprot)=
 600      &         zvdertot(kk,ib,iprot)+
 601      &         zvepstot(icant(iti,j),ib,iprot)
 602             enddo
 603           enddo
 604           do i=1,npair_sc
 605             kk=kk+1
 606             iti=ityp_psc(1,i)
 607             itj=ityp_psc(2,i)
 608             zvdertot(kk,ib,iprot)=
 609      &        zvepstot(icant(iti,itj),ib,iprot)
 610           enddo
 611           do i=1,ntor_var
 612             kk=kk+1
 613             zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
 614           enddo
 615           do i=1,nang_var
 616             kk=kk+1
 617             zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot)
 618           enddo
 619         enddo
 620       enddo
 621 c Molecular properties
 622       do iprot=1,nprot
 623         do inat=1,natlike(iprot)
 624           do ib=1,nbeta(inat+2,iprot)
 625             call propderiv(ib,inat,iprot)
 626             do k=1,natdim(inat,iprot)
 627               do kk=1,ii
 628                 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
 629               enddo
 630               kk=ii
 631               do i=1,nsingle_sc
 632                 kk=kk+1
 633                 iti=ityp_ssc(i)
 634                 rder_all(kk,k,ib,inat,iprot)=
 635      &           reps(icant(iti,iti),k)
 636                 do j=1,ntyp
 637                   rder_all(kk,k,ib,inat,iprot)=
 638      &             rder_all(kk,k,ib,inat,iprot)+
 639      &             reps(icant(iti,j),k)
 640                 enddo
 641               enddo
 642               do i=1,npair_sc
 643                 kk=kk+1
 644                 iti=ityp_psc(1,i)
 645                 itj=ityp_psc(2,i)
 646                 rder_all(kk,k,ib,inat,iprot)=
 647      &           reps(icant(iti,itj),k)
 648               enddo
 649               kk=ii
 650               do i=1,ntor_var
 651                 kk=kk+1
 652                 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
 653               enddo
 654               do i=1,nang_var
 655                 kk=kk+1
 656                 rder_all(kk,k,ib,inat,iprot)=rang(i,k)
 657               enddo
 658             enddo
 659           enddo
 660         enddo
 661       enddo
 662       write (iout,*) "CHECKGRAD: Before entering num targetgfun deriv"
 663       call flush(iout)
 664 c Calculate numerical derivatives and compare with analytical derivatives
 665 c      do i=1,kk
 666       do i=1,n
 667         xi=x(i)
 668         x(i)=xi+xi*delta
 669 c        write (iout,*) (x(k),k=1,kk)
 670         call func1(n,1,x)
 671 c        write (iout,*) "after func1"
 672 c        call flush(iout)
 673 #ifdef NEWCORR
 674 c AL 2/10/2018 Add PMF gradient
 675 c       if(mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 676         if (torsion_pmf.or.turn_pmf.or.eello_pmf.or.angle_PMF)
 677      &    fpotplus = rdif(n,x,gfpot,1)
 678         if (pdbpmf) fpdbplus = maxlik_pdb(n,x,gfpdb,1)
 679 c          write (iout,*) "after rdif"
 680 c          call flush(iout)
 681 c        endif
 682 #endif
 683         write (iout,*) "Variable",i
 684 c Maximum likelihood
 685         do iprot=1,nprot
 686           write (iout,*) "Derivatives of maximum likelihood"
 687           do ib=1,nbeta(1,iprot)
 688 c            write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
 689             pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta)
 690             write (iout,'(2i5,3e14.5)') iprot,ib,
 691      &        sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
 692      &        /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
 693 c            write (iout,*) sumlik(ib,iprot),
 694 c     &        sumlik_orig(ib,iprot)
 695           enddo
 696         enddo
 697 c Heat capacity
 698         do iprot=1,nprot
 699           write (iout,*) "Derivatives of heat capacity"
 700           do ib=1,nbeta(2,iprot)
 701             pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta)
 702             write (iout,'(2i5,3e14.5)') iprot,ib,
 703      &        zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
 704      &        /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
 705 c            write (iout,*) zvtot(ib,ibatch,iprot),
 706 c     &        zvtot_orig(ib,ibatch,iprot)
 707           enddo
 708         enddo
 709 c Numerical derivatives of natlike properties
 710         write (iout,*) "Derivatives of molecular properties"
 711         do iprot=1,nprot
 712           do inat=1,natlike(iprot)
 713             do ib=1,nbeta(inat+2,iprot)
 714               do k=1,natdim(inat,iprot)
 715                 rder_num=(nuave(k,ib,inat,iprot)
 716      &                   -nuave_orig(k,ib,inat,iprot))/(xi*delta)
 717                 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
 718      &             rder_all(kk,k,ib,inat,iprot),
 719      &             rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
 720      &             (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
 721               enddo
 722             enddo
 723           enddo
 724         enddo
 725 #ifdef NEWCORR
 726 c AL 2/10/2018 Add PMF gradient
 727       if (mod_fourier(nloctyp).gt.0.or.mod_fouriertor(nloctyp).gt.0)then
 728           write(iout,*)"Derivative of fpot",(-fpotplus+fpot)/(xi*delta),
 729      &      gfpot(i),(gfpot(i)-(-fpotplus+fpot)/(xi*delta))/gfpot(i)
 730           write(iout,*)"Derivative of fpdb",(fpdbplus-fpdb)/(xi*delta),
 731      &      gfpdb(i),(gfpdb(i)-(fpdbplus-fpdb)/(xi*delta))/gfpdb(i)
 732         endif
 733 #endif
 734         x(i)=xi
 735       enddo
 736       write (iout,*)
 737       do i=1,n
 738         xi=x(i)
 739         x(i)=xi+xi*delta
 740 c        write (iout,*) "i",i," xplus ",x(i),xi
 741         call targetfun(n,x,nf,fplus,uiparm,fdum)
 742         x(i)=xi-xi*delta
 743 c        write (iout,*) "i",i," xminus",x(i),xi
 744         call targetfun(n,x,nf,fminus,uiparm,fdum)
 745         x(i)=xi
 746 c        write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus,
 747 c     &   "obfplus",obfplus," obfminus",obfminus
 748         gnum(i)=0.5d0*(fplus-fminus)/(xi*delta)
 749       enddo
 750       write (iout,*) "Comparison of analytical and numerical gradient"
 751       do i=1,n
 752         write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
 753      &    1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
 754       enddo
 755       write (iout,*)
 756       return
 757       end
 758 c-----------------------------------------------------------------------
 759       subroutine prepare_sc
 760       implicit none
 761       include "DIMENSIONS"
 762       include "DIMENSIONS.ZSCOPT"
 763       include "COMMON.NAMES"
 764       include "COMMON.WEIGHTS"
 765       include "COMMON.INTERACT"
 766       include "COMMON.ENERGIES"
 767       include "COMMON.FFIELD"
 768       include "COMMON.TORSION"
 769       include "COMMON.IOUNITS"
 770       logical lprint
 771       integer i,j,ii,nvarr,iti,itj
 772       double precision rri
 773       double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
 774      &  ratsig1,ratsig2,rsum_max,dwa16,r_augm
 775
 776       lprint=.false.
 777
 778         if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
 779      &    potname(ipot),', exponents are ',expon,2*expon
 780 C-----------------------------------------------------------------------
 781 C Calculate the "working" parameters of SC interactions.
 782         dwa16=2.0d0**(1.0d0/6.0d0)
 783         do i=1,ntyp
 784           do j=i,ntyp
 785             sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 786             sigma(j,i)=sigma(i,j)
 787             rs0(i,j)=dwa16*sigma(i,j)
 788             rs0(j,i)=rs0(i,j)
 789           enddo
 790         enddo
 791         if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 792      &   'Working parameters of the SC interactions:',
 793      &   '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 794      &   '  chi1   ','   chi2   '
 795         do i=1,ntyp
 796           do j=i,ntyp
 797             epsij=eps(i,j)
 798             if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 799               rrij=sigma(i,j)
 800             else
 801               rrij=rr0(i)+rr0(j)
 802             endif
 803             r0(i,j)=rrij
 804             r0(j,i)=rrij
 805             rrij=rrij**expon
 806             epsij=eps(i,j)
 807             sigeps=dsign(1.0D0,epsij)
 808             epsij=dabs(epsij)
 809             aa(i,j)=epsij*rrij*rrij
 810             bb(i,j)=-sigeps*epsij*rrij
 811             aa(j,i)=aa(i,j)
 812             bb(j,i)=bb(i,j)
 813             if (ipot.gt.2) then
 814               sigt1sq=sigma0(i)**2
 815               sigt2sq=sigma0(j)**2
 816               sigii1=sigii(i)
 817               sigii2=sigii(j)
 818               ratsig1=sigt2sq/sigt1sq
 819               ratsig2=1.0D0/ratsig1
 820               chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 821               if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 822               rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 823             else
 824               rsum_max=sigma(i,j)
 825             endif
 826             sigmaii(i,j)=rsum_max
 827             sigmaii(j,i)=rsum_max
 828             if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
 829      &      then
 830               r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 831               augm(i,j)=epsij*r_augm**(2*expon)
 832 c             augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 833               augm(j,i)=augm(i,j)
 834             else
 835               augm(i,j)=0.0D0
 836               augm(j,i)=0.0D0
 837             endif
 838             if (lprint) then
 839               write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 840      &        restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 841      &        sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 842             endif
 843           enddo
 844         enddo
 845       return
 846       end