source/maxlik/src_MD_T_maxlik-NEWCORR/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-5/,
 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       integer ind_shield /25/
 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 .and. k.ne.ind_shield) 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           do i=1,nang_var
 562             kk=kk+1
 563             sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot)
 564           enddo
 565         enddo
 566 c Heat capacity terms
 567         do ib=1,nbeta(2,iprot)
 568           kk=ii
 569           do i=1,nsingle_sc
 570             kk=kk+1
 571             iti=ityp_ssc(i)
 572             zvdertot(kk,ib,iprot)=
 573      &       zvepstot(icant(iti,iti),ib,iprot)
 574             do j=1,ntyp
 575               zvdertot(kk,ib,iprot)=
 576      &         zvdertot(kk,ib,iprot)+
 577      &         zvepstot(icant(iti,j),ib,iprot)
 578             enddo
 579           enddo
 580           do i=1,npair_sc
 581             kk=kk+1
 582             iti=ityp_psc(1,i)
 583             itj=ityp_psc(2,i)
 584             zvdertot(kk,ib,iprot)=
 585      &        zvepstot(icant(iti,itj),ib,iprot)
 586           enddo
 587           do i=1,ntor_var
 588             kk=kk+1
 589             zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
 590           enddo
 591           do i=1,nang_var
 592             kk=kk+1
 593             zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot)
 594           enddo
 595         enddo
 596       enddo
 597 c Molecular properties
 598       do iprot=1,nprot
 599         do inat=1,natlike(iprot)
 600           do ib=1,nbeta(inat+2,iprot)
 601             call propderiv(ib,inat,iprot)
 602             do k=1,natdim(inat,iprot)
 603               do kk=1,ii
 604                 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
 605               enddo
 606               kk=ii
 607               do i=1,nsingle_sc
 608                 kk=kk+1
 609                 iti=ityp_ssc(i)
 610                 rder_all(kk,k,ib,inat,iprot)=
 611      &           reps(icant(iti,iti),k)
 612                 do j=1,ntyp
 613                   rder_all(kk,k,ib,inat,iprot)=
 614      &             rder_all(kk,k,ib,inat,iprot)+
 615      &             reps(icant(iti,j),k)
 616                 enddo
 617               enddo
 618               do i=1,npair_sc
 619                 kk=kk+1
 620                 iti=ityp_psc(1,i)
 621                 itj=ityp_psc(2,i)
 622                 rder_all(kk,k,ib,inat,iprot)=
 623      &           reps(icant(iti,itj),k)
 624               enddo
 625               kk=ii
 626               do i=1,ntor_var
 627                 kk=kk+1
 628                 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
 629               enddo
 630               do i=1,nang_var
 631                 kk=kk+1
 632                 rder_all(kk,k,ib,inat,iprot)=rang(i,k)
 633               enddo
 634             enddo
 635           enddo
 636         enddo
 637       enddo
 638 c Calculate numerical derivatives and compare with analytical derivatives
 639       do i=1,kk
 640         xi=x(i)
 641         x(i)=xi+xi*delta
 642 c        write (iout,*) (x(k),k=1,kk)
 643         call func1(n,1,x)
 644         write (iout,*) "Variable",i
 645 c Maximum likelihood
 646         do iprot=1,nprot
 647           write (iout,*) "Derivatives of maximum likelihood"
 648           do ib=1,nbeta(1,iprot)
 649 c            write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
 650             pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta)
 651             write (iout,'(2i5,3e14.5)') iprot,ib,
 652      &        sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
 653      &        /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
 654 c            write (iout,*) sumlik(ib,iprot),
 655 c     &        sumlik_orig(ib,iprot)
 656           enddo
 657         enddo
 658 c Heat capacity
 659         do iprot=1,nprot
 660           write (iout,*) "Derivatives of heat capacity"
 661           do ib=1,nbeta(2,iprot)
 662             pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta)
 663             write (iout,'(2i5,3e14.5)') iprot,ib,
 664      &        zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
 665      &        /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
 666 c            write (iout,*) zvtot(ib,ibatch,iprot),
 667 c     &        zvtot_orig(ib,ibatch,iprot)
 668           enddo
 669         enddo
 670 c Numerical derivatives of natlike properties
 671         write (iout,*) "Derivatives of molecular properties"
 672         do iprot=1,nprot
 673           do inat=1,natlike(iprot)
 674             do ib=1,nbeta(inat+2,iprot)
 675               do k=1,natdim(inat,iprot)
 676                 rder_num=(nuave(k,ib,inat,iprot)
 677      &                   -nuave_orig(k,ib,inat,iprot))/(xi*delta)
 678                 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
 679      &             rder_all(kk,k,ib,inat,iprot),
 680      &             rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
 681      &             (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
 682               enddo
 683             enddo
 684           enddo
 685         enddo
 686         x(i)=xi
 687       enddo
 688       write (iout,*)
 689       do i=1,n
 690         xi=x(i)
 691         x(i)=xi+xi*delta
 692 c        write (iout,*) "i",i," xplus ",x(i),xi
 693         call targetfun(n,x,nf,fplus,uiparm,fdum)
 694         x(i)=xi-xi*delta
 695 c        write (iout,*) "i",i," xminus",x(i),xi
 696         call targetfun(n,x,nf,fminus,uiparm,fdum)
 697         x(i)=xi
 698 c        write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus,
 699 c     &   "obfplus",obfplus," obfminus",obfminus
 700         gnum(i)=0.5d0*(fplus-fminus)/(xi*delta)
 701       enddo
 702       write (iout,*) "Comparison of analytical and numerical gradient"
 703       do i=1,n
 704         write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
 705      &    1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
 706       enddo
 707       write (iout,*)
 708       return
 709       end
 710 c-----------------------------------------------------------------------
 711       subroutine prepare_sc
 712       implicit none
 713       include "DIMENSIONS"
 714       include "DIMENSIONS.ZSCOPT"
 715       include "COMMON.NAMES"
 716       include "COMMON.WEIGHTS"
 717       include "COMMON.INTERACT"
 718       include "COMMON.ENERGIES"
 719       include "COMMON.FFIELD"
 720       include "COMMON.TORSION"
 721       include "COMMON.IOUNITS"
 722       logical lprint
 723       integer i,j,ii,nvarr,iti,itj
 724       double precision rri
 725       double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
 726      &  ratsig1,ratsig2,rsum_max,dwa16,r_augm
 727
 728       lprint=.false.
 729
 730         if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
 731      &    potname(ipot),', exponents are ',expon,2*expon
 732 C-----------------------------------------------------------------------
 733 C Calculate the "working" parameters of SC interactions.
 734         dwa16=2.0d0**(1.0d0/6.0d0)
 735         do i=1,ntyp
 736           do j=i,ntyp
 737             sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 738             sigma(j,i)=sigma(i,j)
 739             rs0(i,j)=dwa16*sigma(i,j)
 740             rs0(j,i)=rs0(i,j)
 741           enddo
 742         enddo
 743         if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 744      &   'Working parameters of the SC interactions:',
 745      &   '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 746      &   '  chi1   ','   chi2   '
 747         do i=1,ntyp
 748           do j=i,ntyp
 749             epsij=eps(i,j)
 750             if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 751               rrij=sigma(i,j)
 752             else
 753               rrij=rr0(i)+rr0(j)
 754             endif
 755             r0(i,j)=rrij
 756             r0(j,i)=rrij
 757             rrij=rrij**expon
 758             epsij=eps(i,j)
 759             sigeps=dsign(1.0D0,epsij)
 760             epsij=dabs(epsij)
 761             aa(i,j)=epsij*rrij*rrij
 762             bb(i,j)=-sigeps*epsij*rrij
 763             aa(j,i)=aa(i,j)
 764             bb(j,i)=bb(i,j)
 765             if (ipot.gt.2) then
 766               sigt1sq=sigma0(i)**2
 767               sigt2sq=sigma0(j)**2
 768               sigii1=sigii(i)
 769               sigii2=sigii(j)
 770               ratsig1=sigt2sq/sigt1sq
 771               ratsig2=1.0D0/ratsig1
 772               chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 773               if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 774               rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 775             else
 776               rsum_max=sigma(i,j)
 777             endif
 778             sigmaii(i,j)=rsum_max
 779             sigmaii(j,i)=rsum_max
 780             if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
 781      &      then
 782               r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 783               augm(i,j)=epsij*r_augm**(2*expon)
 784 c             augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 785               augm(j,i)=augm(i,j)
 786             else
 787               augm(i,j)=0.0D0
 788               augm(j,i)=0.0D0
 789             endif
 790             if (lprint) then
 791               write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 792      &        restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 793      &        sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 794             endif
 795           enddo
 796         enddo
 797       return
 798       end