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