source/unres/src-HCD-5D/minimize_p.F

   1       subroutine minimize(etot,x,iretcode,nfun)
   2 #ifdef LBFGS
   3       use minima
   4       use inform
   5       use output
   6       use iounit
   7       use scales
   8 #endif
   9       implicit none
  10       include 'DIMENSIONS'
  11 #ifndef LBFGS
  12       integer liv,lv
  13       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
  14 #endif
  15 *********************************************************************
  16 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
  17 * the calling subprogram.                                           *
  18 * when d(i)=1.0, then v(35) is the length of the initial step,      *
  19 * calculated in the usual pythagorean way.                          *
  20 * absolute convergence occurs when the function is within v(31) of  *
  21 * zero. unless you know the minimum value in advance, abs convg     *
  22 * is probably not useful.                                           *
  23 * relative convergence is when the model predicts that the function *
  24 * will decrease by less than v(32)*abs(fun).                        *
  25 *********************************************************************
  26       include 'COMMON.IOUNITS'
  27       include 'COMMON.VAR'
  28       include 'COMMON.GEO'
  29       include 'COMMON.MINIM'
  30       integer icall
  31       common /srutu/ icall
  32 #ifdef LBFGS
  33       double precision grdmin
  34       external funcgrad
  35       external optsave
  36 #else
  37       dimension iv(liv)
  38       double precision v(1:lv)
  39       common /przechowalnia/ v
  40       integer idum
  41       double precision rdum
  42       double precision fdum
  43       external func,gradient,fdum
  44       external func_restr,grad_restr
  45       logical not_done,change,reduce
  46 #endif
  47       double precision x(maxvar),d(maxvar),xx(maxvar)
  48       double precision energia(0:n_ene)
  49       integer i,nvar_restr,nfun,iretcode
  50       double precision etot
  51 c      common /przechowalnia/ v
  52
  53 #ifdef LBFGS
  54       maxiter=maxmin
  55       coordtype='RIGIDBODY'
  56       grdmin=tolf
  57       jout=iout
  58       jprint=print_min_stat
  59       iwrite=0
  60       if (.not. allocated(scale))  allocate (scale(nvar))
  61 c
  62 c     set scaling parameter for function and derivative values;
  63 c     use square root of median eigenvalue of typical Hessian
  64 c
  65       set_scale = .true.
  66 c      nvar = 0
  67       do i = 1, nvar
  68 c         if (use(i)) then
  69 c            do j = 1, 3
  70 c               nvar = nvar + 1
  71                scale(i) = 12.0d0
  72 c            end do
  73 c         end if
  74       end do
  75 c      write (iout,*) "Calling lbfgs"
  76       write (iout,*) 'Calling LBFGS, minimization in angles'
  77       call var_to_geom(nvar,x)
  78       call chainbuild_extconf
  79       call etotal(energia(0))
  80       call enerprint(energia(0))
  81       call lbfgs (nvar,x,etot,grdmin,funcgrad,optsave)
  82       deallocate(scale)
  83       write (iout,*) "Minimized energy",etot
  84 #else
  85       icall = 1
  86
  87       NOT_DONE=.TRUE.
  88
  89 c     DO WHILE (NOT_DONE)
  90
  91       call deflt(2,iv,liv,lv,v)
  92 * 12 means fresh start, dont call deflt
  93       iv(1)=12
  94 * max num of fun calls
  95       if (maxfun.eq.0) maxfun=500
  96       iv(17)=maxfun
  97 * max num of iterations
  98       if (maxmin.eq.0) maxmin=1000
  99       iv(18)=maxmin
 100 * controls output
 101       iv(19)=2
 102 * selects output unit
 103       iv(21)=0
 104       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 105 * 1 means to print out result
 106       iv(22)=print_min_res
 107 * 1 means to print out summary stats
 108       iv(23)=print_min_stat
 109 * 1 means to print initial x and d
 110       iv(24)=print_min_ini
 111 * min val for v(radfac) default is 0.1
 112       v(24)=0.1D0
 113 * max val for v(radfac) default is 4.0
 114       v(25)=2.0D0
 115 c     v(25)=4.0D0
 116 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 117 * the sumsl default is 0.1
 118       v(26)=0.1D0
 119 * false conv if (act fnctn decrease) .lt. v(34)
 120 * the sumsl default is 100*machep
 121       v(34)=v(34)/100.0D0
 122 * absolute convergence
 123       if (tolf.eq.0.0D0) tolf=1.0D-4
 124       v(31)=tolf
 125 * relative convergence
 126       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 127       v(32)=rtolf
 128 * controls initial step size
 129       v(35)=1.0D-1
 130 * large vals of d correspond to small components of step
 131       do i=1,nphi
 132         d(i)=1.0D-1
 133       enddo
 134       do i=nphi+1,nvar
 135         d(i)=1.0D-1
 136       enddo
 137       write (iout,*) 'Calling SUMSL'
 138       call var_to_geom(nvar,x)
 139       call chainbuild_extconf
 140       call intout
 141       call etotal(energia(0))
 142       call enerprint(energia(0))
 143 c     etot = energia(0)
 144       IF (mask_r) THEN
 145        call x2xx(x,xx,nvar_restr)
 146        call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
 147      &                    iv,liv,lv,v,idum,rdum,fdum)
 148        call xx2x(x,xx)
 149       ELSE
 150        call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
 151       ENDIF
 152       etot=v(10)
 153       iretcode=iv(1)
 154 cd    print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
 155 cd    write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
 156 c     call intout
 157 c     change=reduce(x)
 158       call var_to_geom(nvar,x)
 159 c     if (change) then
 160 c       write (iout,'(a)') 'Reduction worked, minimizing again...'
 161 c     else
 162 c       not_done=.false.
 163 c     endif
 164       call chainbuild_extconf
 165 c     call etotal(energia(0))
 166 c     etot=energia(0)
 167 c     call enerprint(energia(0))
 168       nfun=iv(6)
 169
 170 c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
 171
 172 c     ENDDO ! NOT_DONE
 173 #endif
 174       return
 175       end
 176 #ifdef MPI
 177 c----------------------------------------------------------------------------
 178       subroutine ergastulum
 179       implicit none
 180       include 'DIMENSIONS'
 181 #ifdef MPI
 182       include "mpif.h"
 183       double precision time00
 184       integer ierr,ierror
 185 #endif
 186       include 'COMMON.SETUP'
 187       include 'COMMON.DERIV'
 188       include 'COMMON.VAR'
 189       include 'COMMON.IOUNITS'
 190       include 'COMMON.FFIELD'
 191       include 'COMMON.INTERACT'
 192       include 'COMMON.MD'
 193 #ifdef FIVEDIAG
 194        include 'COMMON.LAGRANGE.5diag'
 195 #else
 196        include 'COMMON.LAGRANGE'
 197 #endif
 198       include 'COMMON.TIME1'
 199       double precision z(maxres6),d_a_tmp(maxres6)
 200       double precision edum(0:n_ene),time_order(0:10)
 201       double precision Gcopy(maxres2,maxres2)
 202       common /przechowalnia/ Gcopy
 203       integer icall /0/
 204       integer i,j,iorder
 205 C Workers wait for variables and NF, and NFL from the boss
 206       iorder=0
 207       do while (iorder.ge.0)
 208 c      write (*,*) 'Processor',fg_rank,' CG group',kolor,
 209 c     & ' receives order from Master'
 210 #ifdef MPI
 211         time00=MPI_Wtime()
 212         call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
 213         time_Bcast=time_Bcast+MPI_Wtime()-time00
 214         if (icall.gt.4 .and. iorder.ge.0)
 215      &   time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
 216 #endif
 217        icall=icall+1
 218 c      write (*,*)
 219 c     & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
 220         if (iorder.eq.0) then
 221           call zerograd
 222           call etotal(edum)
 223 c          write (2,*) "After etotal"
 224 c          write (2,*) "dimen",dimen," dimen3",dimen3
 225 c          call flush(2)
 226         else if (iorder.eq.2) then
 227           call zerograd
 228           call etotal_short(edum)
 229 c          write (2,*) "After etotal_short"
 230 c          write (2,*) "dimen",dimen," dimen3",dimen3
 231 c          call flush(2)
 232         else if (iorder.eq.3) then
 233           call zerograd
 234           call etotal_long(edum)
 235 c          write (2,*) "After etotal_long"
 236 c          write (2,*) "dimen",dimen," dimen3",dimen3
 237 c          call flush(2)
 238         else if (iorder.eq.1) then
 239           call sum_gradient
 240 c          write (2,*) "After sum_gradient"
 241 c          write (2,*) "dimen",dimen," dimen3",dimen3
 242 c          call flush(2
 243 #ifndef FIVEDIAG
 244         else if (iorder.eq.4) then
 245           call ginv_mult(z,d_a_tmp)
 246         else if (iorder.eq.5) then
 247 c Setup MD things for a slave
 248           dimen=(nct-nnt+1)+nside
 249           dimen1=(nct-nnt)+(nct-nnt+1)
 250           dimen3=dimen*3
 251 c          write (2,*) "dimen",dimen," dimen3",dimen3
 252 c          call flush(2)
 253           call int_bounds(dimen,igmult_start,igmult_end)
 254           igmult_start=igmult_start-1
 255           call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 256      &     ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 257            my_ng_count=igmult_end-igmult_start
 258           call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 259      &     MPI_INTEGER,FG_COMM,IERROR)
 260 c          write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 261 c          write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 262           myginv_ng_count=maxres2*my_ng_count
 263 c          write (2,*) "igmult_start",igmult_start," igmult_end",
 264 c     &     igmult_end," my_ng_count",my_ng_count
 265 c          call flush(2)
 266           call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 267      &     nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 268           call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 269      &     nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 270 c          write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 271 c          write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 272 c          call flush(2)
 273 c          call MPI_Barrier(FG_COMM,IERROR)
 274           time00=MPI_Wtime()
 275           call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 276      &     nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 277      &     myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 278 #ifdef TIMING
 279           time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 280 #endif
 281           do i=1,dimen
 282             do j=1,2*my_ng_count
 283               ginv(j,i)=gcopy(i,j)
 284             enddo
 285           enddo
 286 c          write (2,*) "dimen",dimen," dimen3",dimen3
 287 c          write (2,*) "End MD setup"
 288 c          call flush(2)
 289 c           write (iout,*) "My chunk of ginv_block"
 290 c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
 291 #endif
 292         else if (iorder.eq.6) then
 293           call int_from_cart1(.false.)
 294         else if (iorder.eq.7) then
 295           call chainbuild_cart
 296         else if (iorder.eq.8) then
 297           call intcartderiv
 298 #ifndef FIVEDIAG
 299         else if (iorder.eq.9) then
 300           call fricmat_mult(z,d_a_tmp)
 301 #endif
 302         else if (iorder.eq.10) then
 303           call setup_fricmat
 304         endif
 305       enddo
 306       write (*,*) 'Processor',fg_rank,' CG group',kolor,
 307      &  ' absolute rank',myrank,' leves ERGASTULUM.'
 308       write(*,*)'Processor',fg_rank,' wait times for respective orders',
 309      & (' order[',i,']',time_order(i),i=0,10)
 310       return
 311       end
 312 #endif
 313 ************************************************************************
 314 #ifdef LBFGS
 315       double precision function funcgrad(x,g)
 316       implicit none
 317       include 'DIMENSIONS'
 318       include 'COMMON.CONTROL'
 319       include 'COMMON.CHAIN'
 320       include 'COMMON.DERIV'
 321       include 'COMMON.VAR'
 322       include 'COMMON.INTERACT'
 323       include 'COMMON.FFIELD'
 324       include 'COMMON.MD'
 325       include 'COMMON.QRESTR'
 326       include 'COMMON.IOUNITS'
 327       include 'COMMON.GEO'
 328       double precision energia(0:n_ene)
 329       double precision x(nvar),g(nvar)
 330       integer i
 331 c     if (jjj.gt.0) then
 332 c      write (iout,*) "in func x"
 333 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 334 c     endif
 335       call var_to_geom(nvar,x)
 336       call zerograd
 337       call chainbuild_extconf
 338       call etotal(energia(0))
 339       call sum_gradient
 340       funcgrad=energia(0)
 341       call cart2intgrad(nvar,g)
 342 C
 343 C Add the components corresponding to local energy terms.
 344 C
 345 c Add the usampl contributions
 346       if (usampl) then
 347          do i=1,nres-3
 348            gloc(i,icg)=gloc(i,icg)+dugamma(i)
 349          enddo
 350          do i=1,nres-2
 351            gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
 352          enddo
 353       endif
 354       do i=1,nvar
 355 cd      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
 356         g(i)=g(i)+gloc(i,icg)
 357       enddo
 358       return
 359       end
 360 #else
 361       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
 362       implicit real*8 (a-h,o-z)
 363       include 'DIMENSIONS'
 364       include 'COMMON.DERIV'
 365       include 'COMMON.IOUNITS'
 366       include 'COMMON.GEO'
 367       common /chuju/ jjj
 368       double precision energia(0:n_ene)
 369       integer jjj
 370       double precision ufparm
 371       external ufparm
 372       integer uiparm(1)
 373       real*8 urparm(1)
 374       dimension x(maxvar)
 375 c     if (jjj.gt.0) then
 376 c      write (iout,*) "in func x"
 377 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 378 c     endif
 379       nfl=nf
 380       icg=mod(nf,2)+1
 381 cd      print *,'func',nf,nfl,icg
 382       call var_to_geom(n,x)
 383       call zerograd
 384       call chainbuild_extconf
 385 cd    write (iout,*) 'ETOTAL called from FUNC'
 386       call etotal(energia(0))
 387       call sum_gradient
 388       f=energia(0)
 389 c     if (jjj.gt.0) then
 390 c      write (iout,*) "upon exit from func"
 391 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 392 c      write (iout,*) 'f=',f
 393 c       jjj=0
 394 c     endif
 395       return
 396       end
 397 ************************************************************************
 398       subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
 399       implicit real*8 (a-h,o-z)
 400       include 'DIMENSIONS'
 401       include 'COMMON.DERIV'
 402       include 'COMMON.IOUNITS'
 403       include 'COMMON.GEO'
 404       common /chuju/ jjj
 405       double precision energia(0:n_ene)
 406       integer jjj
 407       double precision ufparm
 408       external ufparm
 409       integer uiparm(1)
 410       real*8 urparm(1)
 411       dimension x(maxvar)
 412 c     if (jjj.gt.0) then
 413 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 414 c     endif
 415       nfl=nf
 416       icg=mod(nf,2)+1
 417       call var_to_geom_restr(n,x)
 418       call zerograd
 419       call chainbuild_extconf
 420 cd    write (iout,*) 'ETOTAL called from FUNC'
 421       call etotal(energia(0))
 422       call sum_gradient
 423       f=energia(0)
 424 c     if (jjj.gt.0) then
 425 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 426 c       write (iout,*) 'f=',etot
 427 c       jjj=0
 428 c     endif
 429       return
 430       end
 431 c-------------------------------------------------------
 432 #endif
 433       subroutine x2xx(x,xx,n)
 434       implicit real*8 (a-h,o-z)
 435       include 'DIMENSIONS'
 436       include 'COMMON.VAR'
 437       include 'COMMON.CHAIN'
 438       include 'COMMON.INTERACT'
 439       include 'COMMON.IOUNITS'
 440       double precision xx(maxvar),x(maxvar)
 441
 442 c      write (iout,*) "nvar",nvar
 443       do i=1,nvar
 444         varall(i)=x(i)
 445       enddo
 446
 447       ig=0
 448       igall=0
 449       do i=4,nres
 450         igall=igall+1
 451         if (mask_phi(i).eq.1) then
 452           ig=ig+1
 453           xx(ig)=x(igall)
 454         endif
 455       enddo
 456
 457       do i=3,nres
 458         igall=igall+1
 459         if (mask_theta(i).eq.1) then
 460           ig=ig+1
 461           xx(ig)=x(igall)
 462         endif
 463       enddo
 464
 465       do ij=1,2
 466       do i=2,nres-1
 467         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 468           igall=igall+1
 469           if (mask_side(i).eq.1) then
 470             ig=ig+1
 471             xx(ig)=x(igall)
 472 c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
 473 c            write (iout,*) "x",x(igall)," xx",xx(ig)
 474           endif
 475         endif
 476       enddo
 477       enddo
 478
 479       n=ig
 480
 481       return
 482       end
 483
 484       subroutine xx2x(x,xx)
 485       implicit real*8 (a-h,o-z)
 486       include 'DIMENSIONS'
 487       include 'COMMON.VAR'
 488       include 'COMMON.CHAIN'
 489       include 'COMMON.INTERACT'
 490       include 'COMMON.IOUNITS'
 491       double precision xx(maxvar),x(maxvar)
 492
 493       do i=1,nvar
 494         x(i)=varall(i)
 495       enddo
 496
 497       ig=0
 498       igall=0
 499       do i=4,nres
 500         igall=igall+1
 501         if (mask_phi(i).eq.1) then
 502           ig=ig+1
 503           x(igall)=xx(ig)
 504         endif
 505       enddo
 506
 507       do i=3,nres
 508         igall=igall+1
 509         if (mask_theta(i).eq.1) then
 510           ig=ig+1
 511           x(igall)=xx(ig)
 512         endif
 513       enddo
 514
 515       do ij=1,2
 516       do i=2,nres-1
 517         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 518           igall=igall+1
 519           if (mask_side(i).eq.1) then
 520             ig=ig+1
 521             x(igall)=xx(ig)
 522 c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
 523 c            write (iout,*) "x",x(igall)," xx",xx(ig)
 524           endif
 525         endif
 526       enddo
 527       enddo
 528
 529       return
 530       end
 531
 532 c----------------------------------------------------------
 533       subroutine minim_dc(etot,iretcode,nfun)
 534 #ifdef LBFGS
 535       use minima
 536       use inform
 537       use output
 538       use iounit
 539       use scales
 540 #endif
 541       implicit real*8 (a-h,o-z)
 542       include 'DIMENSIONS'
 543 #ifndef LBFGS
 544       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 545 #endif
 546 #ifdef MPI
 547       include 'mpif.h'
 548 #endif
 549       include 'COMMON.CONTROL'
 550       include 'COMMON.SETUP'
 551       include 'COMMON.IOUNITS'
 552       include 'COMMON.VAR'
 553       include 'COMMON.GEO'
 554       include 'COMMON.MINIM'
 555       include 'COMMON.CHAIN'
 556       double precision minval,x(maxvar),d(maxvar),xx(maxvar)
 557 #ifdef LBFGS
 558       double precision grdmin
 559       double precision funcgrad_dc
 560       external funcgrad_dc,optsave
 561 #else
 562       dimension iv(liv)
 563       double precision v(1:lv)
 564       common /przechowalnia/ v
 565       external func_dc,grad_dc,fdum
 566 #endif
 567       double precision g(maxvar),f1
 568       integer nvarx
 569 #ifdef LBFGS
 570       maxiter=maxmin
 571       coordtype='CARTESIAN'
 572       grdmin=tolf
 573       jout=iout
 574       jprint=print_min_stat
 575       iwrite=0
 576 #else
 577       call deflt(2,iv,liv,lv,v)
 578 * 12 means fresh start, dont call deflt
 579       iv(1)=12
 580 * max num of fun calls
 581       if (maxfun.eq.0) maxfun=500
 582       iv(17)=maxfun
 583 * max num of iterations
 584       if (maxmin.eq.0) maxmin=1000
 585       iv(18)=maxmin
 586 * controls output
 587       iv(19)=2
 588 * selects output unit
 589       iv(21)=0
 590       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 591 * 1 means to print out result
 592       iv(22)=print_min_res
 593 * 1 means to print out summary stats
 594       iv(23)=print_min_stat
 595 * 1 means to print initial x and d
 596       iv(24)=print_min_ini
 597 * min val for v(radfac) default is 0.1
 598       v(24)=0.1D0
 599 * max val for v(radfac) default is 4.0
 600       v(25)=2.0D0
 601 c     v(25)=4.0D0
 602 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 603 * the sumsl default is 0.1
 604       v(26)=0.1D0
 605 * false conv if (act fnctn decrease) .lt. v(34)
 606 * the sumsl default is 100*machep
 607       v(34)=v(34)/100.0D0
 608 * absolute convergence
 609       if (tolf.eq.0.0D0) tolf=1.0D-4
 610       v(31)=tolf
 611 * relative convergence
 612       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 613       v(32)=rtolf
 614 * controls initial step size
 615        v(35)=1.0D-1
 616 * large vals of d correspond to small components of step
 617       do i=1,6*nres
 618         d(i)=1.0D-1
 619       enddo
 620 #endif
 621       k=0
 622       do i=1,nres-1
 623         do j=1,3
 624           k=k+1
 625           x(k)=dc(j,i)
 626         enddo
 627       enddo
 628       do i=2,nres-1
 629         if (ialph(i,1).gt.0) then
 630         do j=1,3
 631           k=k+1
 632           x(k)=dc(j,i+nres)
 633         enddo
 634         endif
 635       enddo
 636       nvarx=k
 637       write (iout,*) "Variables set up nvarx",nvarx
 638 #ifdef LBFGS
 639 c
 640 c     From tinker
 641 c
 642 c     perform dynamic allocation of some global arrays
 643 c
 644       if (.not. allocated(scale))  allocate (scale(nvarx))
 645 c
 646 c     set scaling parameter for function and derivative values;
 647 c     use square root of median eigenvalue of typical Hessian
 648 c
 649       set_scale = .true.
 650 c      nvar = 0
 651       do i = 1, nvarx
 652 c         if (use(i)) then
 653 c            do j = 1, 3
 654 c               nvar = nvar + 1
 655                scale(i) = 12.0d0
 656 c            end do
 657 c         end if
 658       end do
 659 c      write (iout,*) "minim_dc Calling lbfgs"
 660       call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
 661       deallocate(scale)
 662 c      write (iout,*) "minim_dc After lbfgs"
 663 #else
 664 c-----
 665 c      write (iout,*) "checkgrad before SUMSL"
 666 c      icheckgrad=1
 667 c      aincr=1.0d-7
 668 c      call exec_checkgrad
 669 c-----
 670
 671       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 672 c-----
 673 c      write (iout,*) "checkgrad after SUMSL"
 674 c      call exec_checkgrad
 675 c-----
 676 #endif
 677       k=0
 678       do i=1,nres-1
 679         do j=1,3
 680           k=k+1
 681           dc(j,i)=x(k)
 682         enddo
 683       enddo
 684       do i=2,nres-1
 685         if (ialph(i,1).gt.0) then
 686         do j=1,3
 687           k=k+1
 688           dc(j,i+nres)=x(k)
 689         enddo
 690         endif
 691       enddo
 692       call chainbuild_cart
 693
 694 cd      call zerograd
 695 cd      nf=0
 696 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 697 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 698 cd
 699 cd      do i=1,k
 700 cd       x(i)=x(i)+1.0D-5
 701 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 702 cd       x(i)=x(i)-1.0D-5
 703 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 704 cd      enddo
 705 #ifndef LBFGS
 706       etot=v(10)
 707       iretcode=iv(1)
 708       nfun=iv(6)
 709 #endif
 710       return
 711       end
 712 #ifdef LBFGS
 713       double precision function funcgrad_dc(x,g)
 714       implicit real*8 (a-h,o-z)
 715       include 'DIMENSIONS'
 716 #ifdef MPI
 717       include 'mpif.h'
 718 #endif
 719       include 'COMMON.SETUP'
 720       include 'COMMON.CHAIN'
 721       include 'COMMON.DERIV'
 722       include 'COMMON.VAR'
 723       include 'COMMON.INTERACT'
 724       include 'COMMON.FFIELD'
 725       include 'COMMON.MD'
 726       include 'COMMON.IOUNITS'
 727       integer k
 728       dimension x(maxvar),g(maxvar)
 729       double precision energia(0:n_ene)
 730       common /gacia/ nf
 731 c
 732       nf=nf+1
 733       k=0
 734       do i=1,nres-1
 735         do j=1,3
 736           k=k+1
 737           dc(j,i)=x(k)
 738         enddo
 739       enddo
 740       do i=2,nres-1
 741         if (ialph(i,1).gt.0) then
 742         do j=1,3
 743           k=k+1
 744           dc(j,i+nres)=x(k)
 745         enddo
 746         endif
 747       enddo
 748       call chainbuild_cart
 749       call zerograd
 750       call etotal(energia(0))
 751 c      write (iout,*) "energia",energia(0)
 752       funcgrad_dc=energia(0)
 753 C
 754 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 755 C
 756       call cartgrad
 757       k=0
 758       do i=1,nres-1
 759         do j=1,3
 760           k=k+1
 761           g(k)=gcart(j,i)
 762         enddo
 763       enddo
 764       do i=2,nres-1
 765         if (ialph(i,1).gt.0) then
 766         do j=1,3
 767           k=k+1
 768           g(k)=gxcart(j,i)
 769         enddo
 770         endif
 771       enddo
 772
 773       return
 774       end
 775 #else
 776       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 777       implicit real*8 (a-h,o-z)
 778       include 'DIMENSIONS'
 779 #ifdef MPI
 780       include 'mpif.h'
 781 #endif
 782       include 'COMMON.SETUP'
 783       include 'COMMON.DERIV'
 784       include 'COMMON.IOUNITS'
 785       include 'COMMON.GEO'
 786       include 'COMMON.CHAIN'
 787       include 'COMMON.VAR'
 788       double precision energia(0:n_ene)
 789       double precision ufparm
 790       external ufparm
 791       integer uiparm(1)
 792       real*8 urparm(1)
 793       dimension x(maxvar)
 794       nfl=nf
 795 cbad      icg=mod(nf,2)+1
 796       icg=1
 797
 798       k=0
 799       do i=1,nres-1
 800         do j=1,3
 801           k=k+1
 802           dc(j,i)=x(k)
 803         enddo
 804       enddo
 805       do i=2,nres-1
 806         if (ialph(i,1).gt.0) then
 807         do j=1,3
 808           k=k+1
 809           dc(j,i+nres)=x(k)
 810         enddo
 811         endif
 812       enddo
 813       call chainbuild_cart
 814
 815       call zerograd
 816       call etotal(energia(0))
 817       f=energia(0)
 818
 819 cd      print *,'func_dc ',nf,nfl,f
 820
 821       return
 822       end
 823
 824       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 825       implicit real*8 (a-h,o-z)
 826       include 'DIMENSIONS'
 827 #ifdef MPI
 828       include 'mpif.h'
 829 #endif
 830       include 'COMMON.SETUP'
 831       include 'COMMON.CHAIN'
 832       include 'COMMON.DERIV'
 833       include 'COMMON.VAR'
 834       include 'COMMON.INTERACT'
 835       include 'COMMON.FFIELD'
 836       include 'COMMON.MD'
 837       include 'COMMON.IOUNITS'
 838       external ufparm
 839       integer uiparm(1),k
 840       double precision urparm(1)
 841       dimension x(maxvar),g(maxvar)
 842 c
 843 c
 844 c
 845 cbad      icg=mod(nf,2)+1
 846       icg=1
 847 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 848       if (nf-nfl+1) 20,30,40
 849    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 850 cd      print *,20
 851       if (nf.eq.0) return
 852       goto 40
 853    30 continue
 854 cd      print *,30
 855       k=0
 856       do i=1,nres-1
 857         do j=1,3
 858           k=k+1
 859           dc(j,i)=x(k)
 860         enddo
 861       enddo
 862       do i=2,nres-1
 863         if (ialph(i,1).gt.0) then
 864         do j=1,3
 865           k=k+1
 866           dc(j,i+nres)=x(k)
 867         enddo
 868         endif
 869       enddo
 870       call chainbuild_cart
 871
 872 C
 873 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 874 C
 875    40 call cartgrad
 876 cd      print *,40
 877       k=0
 878       do i=1,nres-1
 879         do j=1,3
 880           k=k+1
 881           g(k)=gcart(j,i)
 882         enddo
 883       enddo
 884       do i=2,nres-1
 885         if (ialph(i,1).gt.0) then
 886         do j=1,3
 887           k=k+1
 888           g(k)=gxcart(j,i)
 889         enddo
 890         endif
 891       enddo
 892
 893       return
 894       end
 895 #endif