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       write (iout,*) "Variables set up nvarx",nvarx
 639       write (iout,*) "Before energy minimization"
 640       call etotal(energia(0))
 641       call enerprint(energia(0))
 642 #ifdef LBFGS
 643 c
 644 c     From tinker
 645 c
 646 c     perform dynamic allocation of some global arrays
 647 c
 648       if (.not. allocated(scale))  allocate (scale(nvarx))
 649 c
 650 c     set scaling parameter for function and derivative values;
 651 c     use square root of median eigenvalue of typical Hessian
 652 c
 653       set_scale = .true.
 654 c      nvar = 0
 655       do i = 1, nvarx
 656 c         if (use(i)) then
 657 c            do j = 1, 3
 658 c               nvar = nvar + 1
 659                scale(i) = 12.0d0
 660 c            end do
 661 c         end if
 662       end do
 663 c      write (iout,*) "minim_dc Calling lbfgs"
 664       call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
 665       deallocate(scale)
 666 c      write (iout,*) "minim_dc After lbfgs"
 667 #else
 668 c-----
 669 c      write (iout,*) "checkgrad before SUMSL"
 670 c      icheckgrad=1
 671 c      aincr=1.0d-7
 672 c      call exec_checkgrad
 673 c-----
 674
 675       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 676 c-----
 677 c      write (iout,*) "checkgrad after SUMSL"
 678 c      call exec_checkgrad
 679 c-----
 680 #endif
 681       k=0
 682       do i=1,nres-1
 683         do j=1,3
 684           k=k+1
 685           dc(j,i)=x(k)
 686         enddo
 687       enddo
 688       do i=2,nres-1
 689         if (ialph(i,1).gt.0) then
 690         do j=1,3
 691           k=k+1
 692           dc(j,i+nres)=x(k)
 693         enddo
 694         endif
 695       enddo
 696       call chainbuild_cart
 697
 698 cd      call zerograd
 699 cd      nf=0
 700 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 701 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 702 cd
 703 cd      do i=1,k
 704 cd       x(i)=x(i)+1.0D-5
 705 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 706 cd       x(i)=x(i)-1.0D-5
 707 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 708 cd      enddo
 709 #ifndef LBFGS
 710       etot=v(10)
 711       iretcode=iv(1)
 712       nfun=iv(6)
 713 #endif
 714       return
 715       end
 716 #ifdef LBFGS
 717       double precision function funcgrad_dc(x,g)
 718       implicit real*8 (a-h,o-z)
 719       include 'DIMENSIONS'
 720 #ifdef MPI
 721       include 'mpif.h'
 722 #endif
 723       include 'COMMON.SETUP'
 724       include 'COMMON.CHAIN'
 725       include 'COMMON.DERIV'
 726       include 'COMMON.VAR'
 727       include 'COMMON.INTERACT'
 728       include 'COMMON.FFIELD'
 729       include 'COMMON.MD'
 730       include 'COMMON.IOUNITS'
 731       integer k
 732       dimension x(maxvar),g(maxvar)
 733       double precision energia(0:n_ene)
 734       common /gacia/ nf
 735 c
 736       nf=nf+1
 737       k=0
 738       do i=1,nres-1
 739         do j=1,3
 740           k=k+1
 741           dc(j,i)=x(k)
 742         enddo
 743       enddo
 744       do i=2,nres-1
 745         if (ialph(i,1).gt.0) then
 746         do j=1,3
 747           k=k+1
 748           dc(j,i+nres)=x(k)
 749         enddo
 750         endif
 751       enddo
 752       call chainbuild_cart
 753       call zerograd
 754       call etotal(energia(0))
 755 c      write (iout,*) "energia",energia(0)
 756       funcgrad_dc=energia(0)
 757 C
 758 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 759 C
 760       call cartgrad
 761       k=0
 762       do i=1,nres-1
 763         do j=1,3
 764           k=k+1
 765           g(k)=gcart(j,i)
 766         enddo
 767       enddo
 768       do i=2,nres-1
 769         if (ialph(i,1).gt.0) then
 770         do j=1,3
 771           k=k+1
 772           g(k)=gxcart(j,i)
 773         enddo
 774         endif
 775       enddo
 776
 777       return
 778       end
 779 #else
 780       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 781       implicit real*8 (a-h,o-z)
 782       include 'DIMENSIONS'
 783 #ifdef MPI
 784       include 'mpif.h'
 785 #endif
 786       include 'COMMON.SETUP'
 787       include 'COMMON.DERIV'
 788       include 'COMMON.IOUNITS'
 789       include 'COMMON.GEO'
 790       include 'COMMON.CHAIN'
 791       include 'COMMON.VAR'
 792       double precision energia(0:n_ene)
 793       double precision ufparm
 794       external ufparm
 795       integer uiparm(1)
 796       real*8 urparm(1)
 797       dimension x(maxvar)
 798       nfl=nf
 799 cbad      icg=mod(nf,2)+1
 800       icg=1
 801
 802       k=0
 803       do i=1,nres-1
 804         do j=1,3
 805           k=k+1
 806           dc(j,i)=x(k)
 807         enddo
 808       enddo
 809       do i=2,nres-1
 810         if (ialph(i,1).gt.0) then
 811         do j=1,3
 812           k=k+1
 813           dc(j,i+nres)=x(k)
 814         enddo
 815         endif
 816       enddo
 817       call chainbuild_cart
 818
 819       call zerograd
 820       call etotal(energia(0))
 821       f=energia(0)
 822
 823 cd      print *,'func_dc ',nf,nfl,f
 824
 825       return
 826       end
 827
 828       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 829       implicit real*8 (a-h,o-z)
 830       include 'DIMENSIONS'
 831 #ifdef MPI
 832       include 'mpif.h'
 833 #endif
 834       include 'COMMON.SETUP'
 835       include 'COMMON.CHAIN'
 836       include 'COMMON.DERIV'
 837       include 'COMMON.VAR'
 838       include 'COMMON.INTERACT'
 839       include 'COMMON.FFIELD'
 840       include 'COMMON.MD'
 841       include 'COMMON.IOUNITS'
 842       external ufparm
 843       integer uiparm(1),k
 844       double precision urparm(1)
 845       dimension x(maxvar),g(maxvar)
 846 c
 847 c
 848 c
 849 cbad      icg=mod(nf,2)+1
 850       icg=1
 851 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 852       if (nf-nfl+1) 20,30,40
 853    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 854 cd      print *,20
 855       if (nf.eq.0) return
 856       goto 40
 857    30 continue
 858 cd      print *,30
 859       k=0
 860       do i=1,nres-1
 861         do j=1,3
 862           k=k+1
 863           dc(j,i)=x(k)
 864         enddo
 865       enddo
 866       do i=2,nres-1
 867         if (ialph(i,1).gt.0) then
 868         do j=1,3
 869           k=k+1
 870           dc(j,i+nres)=x(k)
 871         enddo
 872         endif
 873       enddo
 874       call chainbuild_cart
 875
 876 C
 877 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 878 C
 879    40 call cartgrad
 880 cd      print *,40
 881       k=0
 882       do i=1,nres-1
 883         do j=1,3
 884           k=k+1
 885           g(k)=gcart(j,i)
 886         enddo
 887       enddo
 888       do i=2,nres-1
 889         if (ialph(i,1).gt.0) then
 890         do j=1,3
 891           k=k+1
 892           g(k)=gxcart(j,i)
 893         enddo
 894         endif
 895       enddo
 896
 897       return
 898       end
 899 #endif