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:11)
 201 c      double precision Gcopy(maxres2,maxres2)
 202 c      common /przechowalnia/ Gcopy
 203       integer icall /0/
 204       integer i,j,iorder,ioverlap(maxres),ioverlap_last
 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         else if (iorder.eq.11) then
 305           call overlap_sc_list(ioverlap,ioverlap_last,.false.)
 306         endif
 307       enddo
 308       write (*,*) 'Processor',fg_rank,' CG group',kolor,
 309      &  ' absolute rank',myrank,' leves ERGASTULUM.'
 310       write(*,*)'Processor',fg_rank,' wait times for respective orders',
 311      & (' order[',i,']',time_order(i),i=0,10)
 312       return
 313       end
 314 #endif
 315 ************************************************************************
 316 #ifdef LBFGS
 317       double precision function funcgrad(x,g)
 318       implicit none
 319       include 'DIMENSIONS'
 320       include 'COMMON.CONTROL'
 321       include 'COMMON.CHAIN'
 322       include 'COMMON.DERIV'
 323       include 'COMMON.VAR'
 324       include 'COMMON.INTERACT'
 325       include 'COMMON.FFIELD'
 326       include 'COMMON.MD'
 327       include 'COMMON.QRESTR'
 328       include 'COMMON.IOUNITS'
 329       include 'COMMON.GEO'
 330       double precision energia(0:n_ene)
 331       double precision x(nvar),g(nvar)
 332       integer i
 333 c     if (jjj.gt.0) then
 334 c      write (iout,*) "in func x"
 335 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 336 c     endif
 337       call var_to_geom(nvar,x)
 338       call zerograd
 339       call chainbuild_extconf
 340       call etotal(energia(0))
 341       call sum_gradient
 342       funcgrad=energia(0)
 343       call cart2intgrad(nvar,g)
 344 C
 345 C Add the components corresponding to local energy terms.
 346 C
 347 c Add the usampl contributions
 348       if (usampl) then
 349          do i=1,nres-3
 350            gloc(i,icg)=gloc(i,icg)+dugamma(i)
 351          enddo
 352          do i=1,nres-2
 353            gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
 354          enddo
 355       endif
 356       do i=1,nvar
 357 cd      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
 358         g(i)=g(i)+gloc(i,icg)
 359       enddo
 360       return
 361       end
 362 #else
 363       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
 364       implicit real*8 (a-h,o-z)
 365       include 'DIMENSIONS'
 366       include 'COMMON.DERIV'
 367       include 'COMMON.IOUNITS'
 368       include 'COMMON.GEO'
 369       common /chuju/ jjj
 370       double precision energia(0:n_ene)
 371       integer jjj
 372       double precision ufparm
 373       external ufparm
 374       integer uiparm(1)
 375       real*8 urparm(1)
 376       dimension x(maxvar)
 377 c     if (jjj.gt.0) then
 378 c      write (iout,*) "in func x"
 379 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 380 c     endif
 381       nfl=nf
 382       icg=mod(nf,2)+1
 383 cd      print *,'func',nf,nfl,icg
 384       call var_to_geom(n,x)
 385       call zerograd
 386       call chainbuild_extconf
 387 cd    write (iout,*) 'ETOTAL called from FUNC'
 388       call etotal(energia(0))
 389       call sum_gradient
 390       f=energia(0)
 391 c     if (jjj.gt.0) then
 392 c      write (iout,*) "upon exit from func"
 393 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 394 c      write (iout,*) 'f=',f
 395 c       jjj=0
 396 c     endif
 397       return
 398       end
 399 ************************************************************************
 400       subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
 401       implicit real*8 (a-h,o-z)
 402       include 'DIMENSIONS'
 403       include 'COMMON.DERIV'
 404       include 'COMMON.IOUNITS'
 405       include 'COMMON.GEO'
 406       common /chuju/ jjj
 407       double precision energia(0:n_ene)
 408       integer jjj
 409       double precision ufparm
 410       external ufparm
 411       integer uiparm(1)
 412       real*8 urparm(1)
 413       dimension x(maxvar)
 414 c     if (jjj.gt.0) then
 415 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 416 c     endif
 417       nfl=nf
 418       icg=mod(nf,2)+1
 419       call var_to_geom_restr(n,x)
 420       call zerograd
 421       call chainbuild_extconf
 422 cd    write (iout,*) 'ETOTAL called from FUNC'
 423       call etotal(energia(0))
 424       call sum_gradient
 425       f=energia(0)
 426 c     if (jjj.gt.0) then
 427 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 428 c       write (iout,*) 'f=',etot
 429 c       jjj=0
 430 c     endif
 431       return
 432       end
 433 c-------------------------------------------------------
 434 #endif
 435       subroutine x2xx(x,xx,n)
 436       implicit real*8 (a-h,o-z)
 437       include 'DIMENSIONS'
 438       include 'COMMON.VAR'
 439       include 'COMMON.CHAIN'
 440       include 'COMMON.INTERACT'
 441       include 'COMMON.IOUNITS'
 442       double precision xx(maxvar),x(maxvar)
 443
 444 c      write (iout,*) "nvar",nvar
 445       do i=1,nvar
 446         varall(i)=x(i)
 447       enddo
 448
 449       ig=0
 450       igall=0
 451       do i=4,nres
 452         igall=igall+1
 453         if (mask_phi(i).eq.1) then
 454           ig=ig+1
 455           xx(ig)=x(igall)
 456         endif
 457       enddo
 458
 459       do i=3,nres
 460         igall=igall+1
 461         if (mask_theta(i).eq.1) then
 462           ig=ig+1
 463           xx(ig)=x(igall)
 464         endif
 465       enddo
 466
 467       do ij=1,2
 468       do i=2,nres-1
 469         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 470           igall=igall+1
 471           if (mask_side(i).eq.1) then
 472             ig=ig+1
 473             xx(ig)=x(igall)
 474 c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
 475 c            write (iout,*) "x",x(igall)," xx",xx(ig)
 476           endif
 477         endif
 478       enddo
 479       enddo
 480
 481       n=ig
 482
 483       return
 484       end
 485
 486       subroutine xx2x(x,xx)
 487       implicit real*8 (a-h,o-z)
 488       include 'DIMENSIONS'
 489       include 'COMMON.VAR'
 490       include 'COMMON.CHAIN'
 491       include 'COMMON.INTERACT'
 492       include 'COMMON.IOUNITS'
 493       double precision xx(maxvar),x(maxvar)
 494
 495       do i=1,nvar
 496         x(i)=varall(i)
 497       enddo
 498
 499       ig=0
 500       igall=0
 501       do i=4,nres
 502         igall=igall+1
 503         if (mask_phi(i).eq.1) then
 504           ig=ig+1
 505           x(igall)=xx(ig)
 506         endif
 507       enddo
 508
 509       do i=3,nres
 510         igall=igall+1
 511         if (mask_theta(i).eq.1) then
 512           ig=ig+1
 513           x(igall)=xx(ig)
 514         endif
 515       enddo
 516
 517       do ij=1,2
 518       do i=2,nres-1
 519         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 520           igall=igall+1
 521           if (mask_side(i).eq.1) then
 522             ig=ig+1
 523             x(igall)=xx(ig)
 524 c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
 525 c            write (iout,*) "x",x(igall)," xx",xx(ig)
 526           endif
 527         endif
 528       enddo
 529       enddo
 530
 531       return
 532       end
 533
 534 c----------------------------------------------------------
 535       subroutine minim_dc(etot,iretcode,nfun)
 536 #ifdef LBFGS
 537       use minima
 538       use inform
 539       use output
 540       use iounit
 541       use scales
 542 #endif
 543       implicit real*8 (a-h,o-z)
 544       include 'DIMENSIONS'
 545 #ifndef LBFGS
 546       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 547 #endif
 548 #ifdef MPI
 549       include 'mpif.h'
 550 #endif
 551       include 'COMMON.CONTROL'
 552       include 'COMMON.SETUP'
 553       include 'COMMON.IOUNITS'
 554       include 'COMMON.VAR'
 555       include 'COMMON.GEO'
 556       include 'COMMON.MINIM'
 557       include 'COMMON.CHAIN'
 558       double precision minval,x(maxvar),d(maxvar),xx(maxvar)
 559 #ifdef LBFGS
 560       double precision grdmin
 561       double precision funcgrad_dc
 562       external funcgrad_dc,optsave
 563 #else
 564       dimension iv(liv)
 565       double precision v(1:lv)
 566       common /przechowalnia/ v
 567       external func_dc,grad_dc,fdum
 568 #endif
 569       double precision g(maxvar),f1
 570       integer nvarx
 571       double precision energia(0:n_ene)
 572 #ifdef LBFGS
 573       maxiter=maxmin
 574       coordtype='CARTESIAN'
 575       grdmin=tolf
 576       jout=iout
 577       jprint=print_min_stat
 578       iwrite=0
 579 #else
 580       call deflt(2,iv,liv,lv,v)
 581 * 12 means fresh start, dont call deflt
 582       iv(1)=12
 583 * max num of fun calls
 584       if (maxfun.eq.0) maxfun=500
 585       iv(17)=maxfun
 586 * max num of iterations
 587       if (maxmin.eq.0) maxmin=1000
 588       iv(18)=maxmin
 589 * controls output
 590       iv(19)=2
 591 * selects output unit
 592       iv(21)=0
 593       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 594 * 1 means to print out result
 595       iv(22)=print_min_res
 596 * 1 means to print out summary stats
 597       iv(23)=print_min_stat
 598 * 1 means to print initial x and d
 599       iv(24)=print_min_ini
 600 * min val for v(radfac) default is 0.1
 601       v(24)=0.1D0
 602 * max val for v(radfac) default is 4.0
 603       v(25)=2.0D0
 604 c     v(25)=4.0D0
 605 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 606 * the sumsl default is 0.1
 607       v(26)=0.1D0
 608 * false conv if (act fnctn decrease) .lt. v(34)
 609 * the sumsl default is 100*machep
 610       v(34)=v(34)/100.0D0
 611 * absolute convergence
 612       if (tolf.eq.0.0D0) tolf=1.0D-4
 613       v(31)=tolf
 614 * relative convergence
 615       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 616       v(32)=rtolf
 617 * controls initial step size
 618        v(35)=1.0D-1
 619 * large vals of d correspond to small components of step
 620       do i=1,6*nres
 621         d(i)=1.0D-1
 622       enddo
 623 #endif
 624       k=0
 625       do i=1,nres-1
 626         do j=1,3
 627           k=k+1
 628           x(k)=dc(j,i)
 629         enddo
 630       enddo
 631       do i=2,nres-1
 632         if (ialph(i,1).gt.0) then
 633         do j=1,3
 634           k=k+1
 635           x(k)=dc(j,i+nres)
 636         enddo
 637         endif
 638       enddo
 639       nvarx=k
 640       call chainbuild_cart
 641       call etotal(energia(0))
 642       call enerprint(energia(0))
 643 #ifdef LBFGS
 644 c
 645 c     From tinker
 646 c
 647 c     perform dynamic allocation of some global arrays
 648 c
 649       if (.not. allocated(scale))  allocate (scale(nvarx))
 650 c
 651 c     set scaling parameter for function and derivative values;
 652 c     use square root of median eigenvalue of typical Hessian
 653 c
 654       set_scale = .true.
 655 c      nvar = 0
 656       do i = 1, nvarx
 657 c         if (use(i)) then
 658 c            do j = 1, 3
 659 c               nvar = nvar + 1
 660                scale(i) = 12.0d0
 661 c            end do
 662 c         end if
 663       end do
 664 c      write (iout,*) "minim_dc Calling lbfgs"
 665       call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
 666       deallocate(scale)
 667 c      write (iout,*) "minim_dc After lbfgs"
 668 #else
 669 c-----
 670 c      write (iout,*) "checkgrad before SUMSL"
 671 c      icheckgrad=1
 672 c      aincr=1.0d-7
 673 c      call exec_checkgrad
 674 c-----
 675
 676       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 677 c-----
 678 c      write (iout,*) "checkgrad after SUMSL"
 679 c      call exec_checkgrad
 680 c-----
 681 #endif
 682       k=0
 683       do i=1,nres-1
 684         do j=1,3
 685           k=k+1
 686           dc(j,i)=x(k)
 687         enddo
 688       enddo
 689       do i=2,nres-1
 690         if (ialph(i,1).gt.0) then
 691         do j=1,3
 692           k=k+1
 693           dc(j,i+nres)=x(k)
 694         enddo
 695         endif
 696       enddo
 697       call chainbuild_cart
 698
 699 cd      call zerograd
 700 cd      nf=0
 701 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 702 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 703 cd
 704 cd      do i=1,k
 705 cd       x(i)=x(i)+1.0D-5
 706 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 707 cd       x(i)=x(i)-1.0D-5
 708 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 709 cd      enddo
 710 #ifndef LBFGS
 711       etot=v(10)
 712       iretcode=iv(1)
 713       nfun=iv(6)
 714 #endif
 715       return
 716       end
 717 #ifdef LBFGS
 718       double precision function funcgrad_dc(x,g)
 719       implicit real*8 (a-h,o-z)
 720       include 'DIMENSIONS'
 721 #ifdef MPI
 722       include 'mpif.h'
 723 #endif
 724       include 'COMMON.SETUP'
 725       include 'COMMON.CHAIN'
 726       include 'COMMON.DERIV'
 727       include 'COMMON.VAR'
 728       include 'COMMON.INTERACT'
 729       include 'COMMON.FFIELD'
 730       include 'COMMON.MD'
 731       include 'COMMON.IOUNITS'
 732       integer k
 733       dimension x(maxvar),g(maxvar)
 734       double precision energia(0:n_ene)
 735       common /gacia/ nf
 736 c
 737       nf=nf+1
 738       k=0
 739       do i=1,nres-1
 740         do j=1,3
 741           k=k+1
 742           dc(j,i)=x(k)
 743         enddo
 744       enddo
 745       do i=2,nres-1
 746         if (ialph(i,1).gt.0) then
 747         do j=1,3
 748           k=k+1
 749           dc(j,i+nres)=x(k)
 750         enddo
 751         endif
 752       enddo
 753       call chainbuild_cart
 754       call zerograd
 755       call etotal(energia(0))
 756 c      write (iout,*) "energia",energia(0)
 757       funcgrad_dc=energia(0)
 758 C
 759 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 760 C
 761       call cartgrad
 762       k=0
 763       do i=1,nres-1
 764         do j=1,3
 765           k=k+1
 766           g(k)=gcart(j,i)
 767         enddo
 768       enddo
 769       do i=2,nres-1
 770         if (ialph(i,1).gt.0) then
 771         do j=1,3
 772           k=k+1
 773           g(k)=gxcart(j,i)
 774         enddo
 775         endif
 776       enddo
 777
 778       return
 779       end
 780 #else
 781       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 782       implicit real*8 (a-h,o-z)
 783       include 'DIMENSIONS'
 784 #ifdef MPI
 785       include 'mpif.h'
 786 #endif
 787       include 'COMMON.SETUP'
 788       include 'COMMON.DERIV'
 789       include 'COMMON.IOUNITS'
 790       include 'COMMON.GEO'
 791       include 'COMMON.CHAIN'
 792       include 'COMMON.VAR'
 793       double precision energia(0:n_ene)
 794       double precision ufparm
 795       external ufparm
 796       integer uiparm(1)
 797       real*8 urparm(1)
 798       dimension x(maxvar)
 799       nfl=nf
 800 cbad      icg=mod(nf,2)+1
 801       icg=1
 802
 803       k=0
 804       do i=1,nres-1
 805         do j=1,3
 806           k=k+1
 807           dc(j,i)=x(k)
 808         enddo
 809       enddo
 810       do i=2,nres-1
 811         if (ialph(i,1).gt.0) then
 812         do j=1,3
 813           k=k+1
 814           dc(j,i+nres)=x(k)
 815         enddo
 816         endif
 817       enddo
 818       call chainbuild_cart
 819
 820       call zerograd
 821       call etotal(energia(0))
 822       f=energia(0)
 823
 824 cd      print *,'func_dc ',nf,nfl,f
 825
 826       return
 827       end
 828
 829       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 830       implicit real*8 (a-h,o-z)
 831       include 'DIMENSIONS'
 832 #ifdef MPI
 833       include 'mpif.h'
 834 #endif
 835       include 'COMMON.SETUP'
 836       include 'COMMON.CHAIN'
 837       include 'COMMON.DERIV'
 838       include 'COMMON.VAR'
 839       include 'COMMON.INTERACT'
 840       include 'COMMON.FFIELD'
 841       include 'COMMON.MD'
 842       include 'COMMON.IOUNITS'
 843       external ufparm
 844       integer uiparm(1),k
 845       double precision urparm(1)
 846       dimension x(maxvar),g(maxvar)
 847 c
 848 c
 849 c
 850 cbad      icg=mod(nf,2)+1
 851       icg=1
 852 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 853       if (nf-nfl+1) 20,30,40
 854    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 855 cd      print *,20
 856       if (nf.eq.0) return
 857       goto 40
 858    30 continue
 859 cd      print *,30
 860       k=0
 861       do i=1,nres-1
 862         do j=1,3
 863           k=k+1
 864           dc(j,i)=x(k)
 865         enddo
 866       enddo
 867       do i=2,nres-1
 868         if (ialph(i,1).gt.0) then
 869         do j=1,3
 870           k=k+1
 871           dc(j,i+nres)=x(k)
 872         enddo
 873         endif
 874       enddo
 875       call chainbuild_cart
 876
 877 C
 878 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 879 C
 880    40 call cartgrad
 881 cd      print *,40
 882       k=0
 883       do i=1,nres-1
 884         do j=1,3
 885           k=k+1
 886           g(k)=gcart(j,i)
 887         enddo
 888       enddo
 889       do i=2,nres-1
 890         if (ialph(i,1).gt.0) then
 891         do j=1,3
 892           k=k+1
 893           g(k)=gxcart(j,i)
 894         enddo
 895         endif
 896       enddo
 897
 898       return
 899       end
 900 #endif