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