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       integer 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       double precision energia(0:n_ene)
 570 #ifdef LBFGS
 571       maxiter=maxmin
 572       coordtype='CARTESIAN'
 573       grdmin=tolf
 574       jout=iout
 575       jprint=print_min_stat
 576       iwrite=0
 577 #else
 578       call deflt(2,iv,liv,lv,v)
 579 * 12 means fresh start, dont call deflt
 580       iv(1)=12
 581 * max num of fun calls
 582       if (maxfun.eq.0) maxfun=500
 583       iv(17)=maxfun
 584 * max num of iterations
 585       if (maxmin.eq.0) maxmin=1000
 586       iv(18)=maxmin
 587 * controls output
 588       iv(19)=2
 589 * selects output unit
 590       iv(21)=0
 591       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 592 * 1 means to print out result
 593       iv(22)=print_min_res
 594 * 1 means to print out summary stats
 595       iv(23)=print_min_stat
 596 * 1 means to print initial x and d
 597       iv(24)=print_min_ini
 598 * min val for v(radfac) default is 0.1
 599       v(24)=0.1D0
 600 * max val for v(radfac) default is 4.0
 601       v(25)=2.0D0
 602 c     v(25)=4.0D0
 603 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 604 * the sumsl default is 0.1
 605       v(26)=0.1D0
 606 * false conv if (act fnctn decrease) .lt. v(34)
 607 * the sumsl default is 100*machep
 608       v(34)=v(34)/100.0D0
 609 * absolute convergence
 610       if (tolf.eq.0.0D0) tolf=1.0D-4
 611       v(31)=tolf
 612 * relative convergence
 613       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 614       v(32)=rtolf
 615 * controls initial step size
 616        v(35)=1.0D-1
 617 * large vals of d correspond to small components of step
 618       do i=1,6*nres
 619         d(i)=1.0D-1
 620       enddo
 621 #endif
 622       k=0
 623       do i=1,nres-1
 624         do j=1,3
 625           k=k+1
 626           x(k)=dc(j,i)
 627         enddo
 628       enddo
 629       do i=2,nres-1
 630         if (ialph(i,1).gt.0) then
 631         do j=1,3
 632           k=k+1
 633           x(k)=dc(j,i+nres)
 634         enddo
 635         endif
 636       enddo
 637       nvarx=k
 638       call chainbuild_cart
 639       call etotal(energia(0))
 640       call enerprint(energia(0))
 641 #ifdef LBFGS
 642 c
 643 c     From tinker
 644 c
 645 c     perform dynamic allocation of some global arrays
 646 c
 647       if (.not. allocated(scale))  allocate (scale(nvarx))
 648 c
 649 c     set scaling parameter for function and derivative values;
 650 c     use square root of median eigenvalue of typical Hessian
 651 c
 652       set_scale = .true.
 653 c      nvar = 0
 654       do i = 1, nvarx
 655 c         if (use(i)) then
 656 c            do j = 1, 3
 657 c               nvar = nvar + 1
 658                scale(i) = 12.0d0
 659 c            end do
 660 c         end if
 661       end do
 662 c      write (iout,*) "minim_dc Calling lbfgs"
 663       call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
 664       deallocate(scale)
 665 c      write (iout,*) "minim_dc After lbfgs"
 666 #else
 667 c-----
 668 c      write (iout,*) "checkgrad before SUMSL"
 669 c      icheckgrad=1
 670 c      aincr=1.0d-7
 671 c      call exec_checkgrad
 672 c-----
 673
 674       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 675 c-----
 676 c      write (iout,*) "checkgrad after SUMSL"
 677 c      call exec_checkgrad
 678 c-----
 679 #endif
 680       k=0
 681       do i=1,nres-1
 682         do j=1,3
 683           k=k+1
 684           dc(j,i)=x(k)
 685         enddo
 686       enddo
 687       do i=2,nres-1
 688         if (ialph(i,1).gt.0) then
 689         do j=1,3
 690           k=k+1
 691           dc(j,i+nres)=x(k)
 692         enddo
 693         endif
 694       enddo
 695       call chainbuild_cart
 696
 697 cd      call zerograd
 698 cd      nf=0
 699 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 700 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 701 cd
 702 cd      do i=1,k
 703 cd       x(i)=x(i)+1.0D-5
 704 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 705 cd       x(i)=x(i)-1.0D-5
 706 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 707 cd      enddo
 708 #ifndef LBFGS
 709       etot=v(10)
 710       iretcode=iv(1)
 711       nfun=iv(6)
 712 #endif
 713       return
 714       end
 715 #ifdef LBFGS
 716       double precision function funcgrad_dc(x,g)
 717       implicit real*8 (a-h,o-z)
 718       include 'DIMENSIONS'
 719 #ifdef MPI
 720       include 'mpif.h'
 721 #endif
 722       include 'COMMON.SETUP'
 723       include 'COMMON.CHAIN'
 724       include 'COMMON.DERIV'
 725       include 'COMMON.VAR'
 726       include 'COMMON.INTERACT'
 727       include 'COMMON.FFIELD'
 728       include 'COMMON.MD'
 729       include 'COMMON.IOUNITS'
 730       integer k
 731       dimension x(maxvar),g(maxvar)
 732       double precision energia(0:n_ene)
 733       common /gacia/ nf
 734 c
 735       nf=nf+1
 736       k=0
 737       do i=1,nres-1
 738         do j=1,3
 739           k=k+1
 740           dc(j,i)=x(k)
 741         enddo
 742       enddo
 743       do i=2,nres-1
 744         if (ialph(i,1).gt.0) then
 745         do j=1,3
 746           k=k+1
 747           dc(j,i+nres)=x(k)
 748         enddo
 749         endif
 750       enddo
 751       call chainbuild_cart
 752       call zerograd
 753       call etotal(energia(0))
 754 c      write (iout,*) "energia",energia(0)
 755       funcgrad_dc=energia(0)
 756 C
 757 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 758 C
 759       call cartgrad
 760       k=0
 761       do i=1,nres-1
 762         do j=1,3
 763           k=k+1
 764           g(k)=gcart(j,i)
 765         enddo
 766       enddo
 767       do i=2,nres-1
 768         if (ialph(i,1).gt.0) then
 769         do j=1,3
 770           k=k+1
 771           g(k)=gxcart(j,i)
 772         enddo
 773         endif
 774       enddo
 775
 776       return
 777       end
 778 #else
 779       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 780       implicit real*8 (a-h,o-z)
 781       include 'DIMENSIONS'
 782 #ifdef MPI
 783       include 'mpif.h'
 784 #endif
 785       include 'COMMON.SETUP'
 786       include 'COMMON.DERIV'
 787       include 'COMMON.IOUNITS'
 788       include 'COMMON.GEO'
 789       include 'COMMON.CHAIN'
 790       include 'COMMON.VAR'
 791       double precision energia(0:n_ene)
 792       double precision ufparm
 793       external ufparm
 794       integer uiparm(1)
 795       real*8 urparm(1)
 796       dimension x(maxvar)
 797       nfl=nf
 798 cbad      icg=mod(nf,2)+1
 799       icg=1
 800
 801       k=0
 802       do i=1,nres-1
 803         do j=1,3
 804           k=k+1
 805           dc(j,i)=x(k)
 806         enddo
 807       enddo
 808       do i=2,nres-1
 809         if (ialph(i,1).gt.0) then
 810         do j=1,3
 811           k=k+1
 812           dc(j,i+nres)=x(k)
 813         enddo
 814         endif
 815       enddo
 816       call chainbuild_cart
 817
 818       call zerograd
 819       call etotal(energia(0))
 820       f=energia(0)
 821
 822 cd      print *,'func_dc ',nf,nfl,f
 823
 824       return
 825       end
 826
 827       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 828       implicit real*8 (a-h,o-z)
 829       include 'DIMENSIONS'
 830 #ifdef MPI
 831       include 'mpif.h'
 832 #endif
 833       include 'COMMON.SETUP'
 834       include 'COMMON.CHAIN'
 835       include 'COMMON.DERIV'
 836       include 'COMMON.VAR'
 837       include 'COMMON.INTERACT'
 838       include 'COMMON.FFIELD'
 839       include 'COMMON.MD'
 840       include 'COMMON.IOUNITS'
 841       external ufparm
 842       integer uiparm(1),k
 843       double precision urparm(1)
 844       dimension x(maxvar),g(maxvar)
 845 c
 846 c
 847 c
 848 cbad      icg=mod(nf,2)+1
 849       icg=1
 850 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 851       if (nf-nfl+1) 20,30,40
 852    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 853 cd      print *,20
 854       if (nf.eq.0) return
 855       goto 40
 856    30 continue
 857 cd      print *,30
 858       k=0
 859       do i=1,nres-1
 860         do j=1,3
 861           k=k+1
 862           dc(j,i)=x(k)
 863         enddo
 864       enddo
 865       do i=2,nres-1
 866         if (ialph(i,1).gt.0) then
 867         do j=1,3
 868           k=k+1
 869           dc(j,i+nres)=x(k)
 870         enddo
 871         endif
 872       enddo
 873       call chainbuild_cart
 874
 875 C
 876 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 877 C
 878    40 call cartgrad
 879 cd      print *,40
 880       k=0
 881       do i=1,nres-1
 882         do j=1,3
 883           k=k+1
 884           g(k)=gcart(j,i)
 885         enddo
 886       enddo
 887       do i=2,nres-1
 888         if (ialph(i,1).gt.0) then
 889         do j=1,3
 890           k=k+1
 891           g(k)=gxcart(j,i)
 892         enddo
 893         endif
 894       enddo
 895
 896       return
 897       end
 898 #endif