minimize_p.F

   1       subroutine minimize(etot,x,iretcode,nfun)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
   5 *********************************************************************
   6 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
   7 * the calling subprogram.                                           *
   8 * when d(i)=1.0, then v(35) is the length of the initial step,      *
   9 * calculated in the usual pythagorean way.                          *
  10 * absolute convergence occurs when the function is within v(31) of  *
  11 * zero. unless you know the minimum value in advance, abs convg     *
  12 * is probably not useful.                                           *
  13 * relative convergence is when the model predicts that the function *
  14 * will decrease by less than v(32)*abs(fun).                        *
  15 *********************************************************************
  16       include 'COMMON.IOUNITS'
  17       include 'COMMON.VAR'
  18       include 'COMMON.GEO'
  19       include 'COMMON.MINIM'
  20       common /srutu/ icall
  21       dimension iv(liv)
  22       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
  23       double precision energia(0:n_ene)
  24       external func,gradient,fdum
  25       external func_restr,grad_restr
  26       logical not_done,change,reduce
  27 c      common /przechowalnia/ v
  28
  29       icall = 1
  30 C Following 3 line for diagnostics; comment out if not needed
  31       write(iout,*) "Enter MINIMIZE liv",liv," lv",lv
  32       write (iout,*) "Coordinates before minimization"
  33       call intout
  34       call deflt(2,iv,liv,lv,v)
  35 * 12 means fresh start, dont call deflt
  36       iv(1)=12
  37 * max num of fun calls
  38       if (maxfun.eq.0) maxfun=500
  39       iv(17)=maxfun
  40 * max num of iterations
  41       if (maxmin.eq.0) maxmin=1000
  42       iv(18)=maxmin
  43 * controls output
  44       iv(19)=2
  45 * selects output unit
  46       iv(21)=0
  47       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
  48 * 1 means to print out result
  49       iv(22)=print_min_res
  50 * 1 means to print out summary stats
  51       iv(23)=print_min_stat
  52 * 1 means to print initial x and d
  53       iv(24)=print_min_ini
  54 * min val for v(radfac) default is 0.1
  55       v(24)=0.1D0
  56 * max val for v(radfac) default is 4.0
  57       v(25)=2.0D0
  58 c     v(25)=4.0D0
  59 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
  60 * the sumsl default is 0.1
  61       v(26)=0.1D0
  62 * false conv if (act fnctn decrease) .lt. v(34)
  63 * the sumsl default is 100*machep
  64       v(34)=v(34)/100.0D0
  65 * absolute convergence
  66       if (tolf.eq.0.0D0) tolf=1.0D-4
  67       v(31)=tolf
  68 * relative convergence
  69       if (rtolf.eq.0.0D0) rtolf=1.0D-4
  70       v(32)=rtolf
  71 * controls initial step size
  72        v(35)=1.0D-1
  73 * large vals of d correspond to small components of step
  74       do i=1,nphi
  75         d(i)=1.0D-1
  76       enddo
  77       do i=nphi+1,nvar
  78         d(i)=1.0D-1
  79       enddo
  80 cd    print *,'Calling SUMSL'
  81       IF (mask_r) THEN
  82        call x2xx(x,xx,nvar_restr)
  83        call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
  84      &                    iv,liv,lv,v,idum,rdum,fdum)
  85        call xx2x(x,xx)
  86       ELSE
  87        call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
  88       ENDIF
  89       etot=v(10)
  90       iretcode=iv(1)
  91 cd    print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
  92       write (iout,'(/a,i4,a,f15.5/)') 'SUMSL return code:',iv(1),
  93      &   " energy",etot
  94       call var_to_geom(nvar,x)
  95 c      call chainbuild
  96 c      call etotal(energia(0))
  97 c      etot=energia(0)
  98 c      call enerprint(energia(0))
  99       nfun=iv(6)
 100
 101       return
 102       end
 103 #ifdef MPI
 104 c----------------------------------------------------------------------------
 105       subroutine ergastulum
 106       implicit real*8 (a-h,o-z)
 107       include 'DIMENSIONS'
 108 #ifdef MPI
 109       include "mpif.h"
 110 #endif
 111       include 'COMMON.SETUP'
 112       include 'COMMON.DERIV'
 113       include 'COMMON.VAR'
 114       include 'COMMON.IOUNITS'
 115       include 'COMMON.FFIELD'
 116       include 'COMMON.INTERACT'
 117       include 'COMMON.MD_'
 118       include 'COMMON.TIME1'
 119       double precision z(maxres6),d_a_tmp(maxres6)
 120       double precision edum(0:n_ene),time_order(0:10)
 121       double precision Gcopy(maxres2,maxres2)
 122       common /przechowalnia/ Gcopy
 123       integer icall /0/
 124 C Workers wait for variables and NF, and NFL from the boss
 125       iorder=0
 126       do while (iorder.ge.0)
 127 c      write (*,*) 'Processor',fg_rank,' CG group',kolor,
 128 c     & ' receives order from Master'
 129         time00=MPI_Wtime()
 130         call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
 131         time_Bcast=time_Bcast+MPI_Wtime()-time00
 132         if (icall.gt.4 .and. iorder.ge.0)
 133      &   time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
 134        icall=icall+1
 135 c      write (*,*)
 136 c     & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
 137         if (iorder.eq.0) then
 138           call zerograd
 139           call etotal(edum)
 140 c          write (2,*) "After etotal"
 141 c          write (2,*) "dimen",dimen," dimen3",dimen3
 142 c          call flush(2)
 143         else if (iorder.eq.2) then
 144           call zerograd
 145 cmd          call etotal_short(edum)
 146 c          write (2,*) "After etotal_short"
 147 c          write (2,*) "dimen",dimen," dimen3",dimen3
 148 c          call flush(2)
 149         else if (iorder.eq.3) then
 150           call zerograd
 151 cmd          call etotal_long(edum)
 152 c          write (2,*) "After etotal_long"
 153 c          write (2,*) "dimen",dimen," dimen3",dimen3
 154 c          call flush(2)
 155         else if (iorder.eq.1) then
 156           call sum_gradient
 157 c          write (2,*) "After sum_gradient"
 158 c          write (2,*) "dimen",dimen," dimen3",dimen3
 159 c          call flush(2)
 160         else if (iorder.eq.4) then
 161 cmd          call ginv_mult(z,d_a_tmp)
 162         else if (iorder.eq.5) then
 163 c Setup MD things for a slave
 164           dimen=(nct-nnt+1)+nside
 165           dimen1=(nct-nnt)+(nct-nnt+1)
 166           dimen3=dimen*3
 167 c          write (2,*) "dimen",dimen," dimen3",dimen3
 168 c          call flush(2)
 169           call int_bounds(dimen,igmult_start,igmult_end)
 170           igmult_start=igmult_start-1
 171           call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 172      &     ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 173            my_ng_count=igmult_end-igmult_start
 174           call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 175      &     MPI_INTEGER,FG_COMM,IERROR)
 176 c          write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 177 c          write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 178           myginv_ng_count=maxres2*my_ng_count
 179 c          write (2,*) "igmult_start",igmult_start," igmult_end",
 180 c     &     igmult_end," my_ng_count",my_ng_count
 181 c          call flush(2)
 182           call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 183      &     nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 184           call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 185      &     nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 186 c          write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 187 c          write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 188 c          call flush(2)
 189 c          call MPI_Barrier(FG_COMM,IERROR)
 190           time00=MPI_Wtime()
 191           call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 192      &     nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 193      &     myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 194 #ifdef TIMING
 195           time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 196 #endif
 197           do i=1,dimen
 198             do j=1,2*my_ng_count
 199               ginv(j,i)=gcopy(i,j)
 200             enddo
 201           enddo
 202 c          write (2,*) "dimen",dimen," dimen3",dimen3
 203 c          write (2,*) "End MD setup"
 204 c          call flush(2)
 205 c           write (iout,*) "My chunk of ginv_block"
 206 c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
 207         else if (iorder.eq.6) then
 208           call int_from_cart1(.false.)
 209         else if (iorder.eq.7) then
 210           call chainbuild_cart
 211         else if (iorder.eq.8) then
 212           call intcartderiv
 213         else if (iorder.eq.9) then
 214 cmd          call fricmat_mult(z,d_a_tmp)
 215         else if (iorder.eq.10) then
 216 cmd          call setup_fricmat
 217         endif
 218       enddo
 219       write (*,*) 'Processor',fg_rank,' CG group',kolor,
 220      &  ' absolute rank',myrank,' leves ERGASTULUM.'
 221       write(*,*)'Processor',fg_rank,' wait times for respective orders',
 222      & (' order[',i,']',time_order(i),i=0,10)
 223       return
 224       end
 225 #endif
 226 ************************************************************************
 227       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
 228       implicit real*8 (a-h,o-z)
 229       include 'DIMENSIONS'
 230       include 'COMMON.DERIV'
 231       include 'COMMON.IOUNITS'
 232       include 'COMMON.GEO'
 233       common /chuju/ jjj
 234       double precision energia(0:n_ene)
 235       integer jjj
 236       double precision ufparm
 237       external ufparm
 238       integer uiparm(1)
 239       real*8 urparm(1)
 240       dimension x(maxvar)
 241 c     if (jjj.gt.0) then
 242 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 243 c     endif
 244       nfl=nf
 245       icg=mod(nf,2)+1
 246 cd      print *,'func',nf,nfl,icg
 247       call var_to_geom(n,x)
 248       call zerograd
 249       call chainbuild
 250 cd    write (iout,*) 'ETOTAL called from FUNC'
 251       call etotal(energia(0))
 252       call sum_gradient
 253       f=energia(0)
 254 c     if (jjj.gt.0) then
 255 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 256 c       write (iout,*) 'f=',etot
 257 c       jjj=0
 258 c     endif
 259       return
 260       end
 261 ************************************************************************
 262       subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
 263       implicit real*8 (a-h,o-z)
 264       include 'DIMENSIONS'
 265       include 'COMMON.DERIV'
 266       include 'COMMON.IOUNITS'
 267       include 'COMMON.GEO'
 268       common /chuju/ jjj
 269       double precision energia(0:n_ene)
 270       integer jjj
 271       double precision ufparm
 272       external ufparm
 273       integer uiparm(1)
 274       real*8 urparm(1)
 275       dimension x(maxvar)
 276 c     if (jjj.gt.0) then
 277 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 278 c     endif
 279       nfl=nf
 280       icg=mod(nf,2)+1
 281       call var_to_geom_restr(n,x)
 282       call zerograd
 283       call chainbuild
 284 cd    write (iout,*) 'ETOTAL called from FUNC'
 285       call etotal(energia(0))
 286       call sum_gradient
 287       f=energia(0)
 288 c     if (jjj.gt.0) then
 289 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 290 c       write (iout,*) 'f=',etot
 291 c       jjj=0
 292 c     endif
 293       return
 294       end
 295 c-------------------------------------------------------
 296       subroutine x2xx(x,xx,n)
 297       implicit real*8 (a-h,o-z)
 298       include 'DIMENSIONS'
 299       include 'COMMON.VAR'
 300       include 'COMMON.CHAIN'
 301       include 'COMMON.INTERACT'
 302       double precision xx(maxvar),x(maxvar)
 303
 304       do i=1,nvar
 305         varall(i)=x(i)
 306       enddo
 307
 308       ig=0
 309       igall=0
 310       do i=4,nres
 311         igall=igall+1
 312         if (mask_phi(i).eq.1) then
 313           ig=ig+1
 314           xx(ig)=x(igall)
 315         endif
 316       enddo
 317
 318       do i=3,nres
 319         igall=igall+1
 320         if (mask_theta(i).eq.1) then
 321           ig=ig+1
 322           xx(ig)=x(igall)
 323         endif
 324       enddo
 325
 326       do ij=1,2
 327       do i=2,nres-1
 328         if (itype(i).ne.10) then
 329           igall=igall+1
 330           if (mask_side(i).eq.1) then
 331             ig=ig+1
 332             xx(ig)=x(igall)
 333           endif
 334         endif
 335       enddo
 336       enddo
 337
 338       n=ig
 339
 340       return
 341       end
 342
 343       subroutine xx2x(x,xx)
 344       implicit real*8 (a-h,o-z)
 345       include 'DIMENSIONS'
 346       include 'COMMON.VAR'
 347       include 'COMMON.CHAIN'
 348       include 'COMMON.INTERACT'
 349       double precision xx(maxvar),x(maxvar)
 350
 351       do i=1,nvar
 352         x(i)=varall(i)
 353       enddo
 354
 355       ig=0
 356       igall=0
 357       do i=4,nres
 358         igall=igall+1
 359         if (mask_phi(i).eq.1) then
 360           ig=ig+1
 361           x(igall)=xx(ig)
 362         endif
 363       enddo
 364
 365       do i=3,nres
 366         igall=igall+1
 367         if (mask_theta(i).eq.1) then
 368           ig=ig+1
 369           x(igall)=xx(ig)
 370         endif
 371       enddo
 372
 373       do ij=1,2
 374       do i=2,nres-1
 375         if (itype(i).ne.10) then
 376           igall=igall+1
 377           if (mask_side(i).eq.1) then
 378             ig=ig+1
 379             x(igall)=xx(ig)
 380           endif
 381         endif
 382       enddo
 383       enddo
 384
 385       return
 386       end
 387
 388 c----------------------------------------------------------
 389       subroutine minim_dc(etot,iretcode,nfun)
 390       implicit real*8 (a-h,o-z)
 391       include 'DIMENSIONS'
 392       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 393 #ifdef MPI
 394       include 'mpif.h'
 395 #endif
 396       include 'COMMON.SETUP'
 397       include 'COMMON.IOUNITS'
 398       include 'COMMON.VAR'
 399       include 'COMMON.GEO'
 400       include 'COMMON.MINIM'
 401       include 'COMMON.CHAIN'
 402       dimension iv(liv)
 403       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
 404 c      common /przechowalnia/ v
 405
 406       double precision energia(0:n_ene)
 407       external func_dc,grad_dc,fdum
 408       logical not_done,change,reduce
 409       double precision g(maxvar),f1
 410
 411       call deflt(2,iv,liv,lv,v)
 412 * 12 means fresh start, dont call deflt
 413       iv(1)=12
 414 * max num of fun calls
 415       if (maxfun.eq.0) maxfun=500
 416       iv(17)=maxfun
 417 * max num of iterations
 418       if (maxmin.eq.0) maxmin=1000
 419       iv(18)=maxmin
 420 * controls output
 421       iv(19)=2
 422 * selects output unit
 423       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 424 * 1 means to print out result
 425       iv(22)=print_min_res
 426 * 1 means to print out summary stats
 427       iv(23)=print_min_stat
 428 * 1 means to print initial x and d
 429       iv(24)=print_min_ini
 430 * min val for v(radfac) default is 0.1
 431       v(24)=0.1D0
 432 * max val for v(radfac) default is 4.0
 433       v(25)=2.0D0
 434 c     v(25)=4.0D0
 435 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 436 * the sumsl default is 0.1
 437       v(26)=0.1D0
 438 * false conv if (act fnctn decrease) .lt. v(34)
 439 * the sumsl default is 100*machep
 440       v(34)=v(34)/100.0D0
 441 * absolute convergence
 442       if (tolf.eq.0.0D0) tolf=1.0D-4
 443       v(31)=tolf
 444 * relative convergence
 445       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 446       v(32)=rtolf
 447 * controls initial step size
 448        v(35)=1.0D-1
 449 * large vals of d correspond to small components of step
 450       do i=1,6*nres
 451         d(i)=1.0D-1
 452       enddo
 453
 454       k=0
 455       do i=1,nres-1
 456         do j=1,3
 457           k=k+1
 458           x(k)=dc(j,i)
 459         enddo
 460       enddo
 461       do i=2,nres-1
 462         if (ialph(i,1).gt.0) then
 463         do j=1,3
 464           k=k+1
 465           x(k)=dc(j,i+nres)
 466         enddo
 467         endif
 468       enddo
 469
 470       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 471
 472       k=0
 473       do i=1,nres-1
 474         do j=1,3
 475           k=k+1
 476           dc(j,i)=x(k)
 477         enddo
 478       enddo
 479       do i=2,nres-1
 480         if (ialph(i,1).gt.0) then
 481         do j=1,3
 482           k=k+1
 483           dc(j,i+nres)=x(k)
 484         enddo
 485         endif
 486       enddo
 487       call chainbuild_cart
 488
 489 cd      call zerograd
 490 cd      nf=0
 491 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 492 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 493 cd
 494 cd      do i=1,k
 495 cd       x(i)=x(i)+1.0D-5
 496 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 497 cd       x(i)=x(i)-1.0D-5
 498 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 499 cd      enddo
 500
 501       etot=v(10)
 502       iretcode=iv(1)
 503       nfun=iv(6)
 504       return
 505       end
 506
 507       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 508       implicit real*8 (a-h,o-z)
 509       include 'DIMENSIONS'
 510 #ifdef MPI
 511       include 'mpif.h'
 512 #endif
 513       include 'COMMON.SETUP'
 514       include 'COMMON.DERIV'
 515       include 'COMMON.IOUNITS'
 516       include 'COMMON.GEO'
 517       include 'COMMON.CHAIN'
 518       include 'COMMON.VAR'
 519       double precision energia(0:n_ene)
 520       double precision ufparm
 521       external ufparm
 522       integer uiparm(1)
 523       real*8 urparm(1)
 524       dimension x(maxvar)
 525       nfl=nf
 526 cbad      icg=mod(nf,2)+1
 527       icg=1
 528
 529       k=0
 530       do i=1,nres-1
 531         do j=1,3
 532           k=k+1
 533           dc(j,i)=x(k)
 534         enddo
 535       enddo
 536       do i=2,nres-1
 537         if (ialph(i,1).gt.0) then
 538         do j=1,3
 539           k=k+1
 540           dc(j,i+nres)=x(k)
 541         enddo
 542         endif
 543       enddo
 544       call chainbuild_cart
 545
 546       call zerograd
 547       call etotal(energia(0))
 548       f=energia(0)
 549
 550 cd      print *,'func_dc ',nf,nfl,f
 551
 552       return
 553       end
 554
 555       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 556       implicit real*8 (a-h,o-z)
 557       include 'DIMENSIONS'
 558 #ifdef MPI
 559       include 'mpif.h'
 560 #endif
 561       include 'COMMON.SETUP'
 562       include 'COMMON.CHAIN'
 563       include 'COMMON.DERIV'
 564       include 'COMMON.VAR'
 565       include 'COMMON.INTERACT'
 566       include 'COMMON.FFIELD'
 567       include 'COMMON.MD_'
 568       include 'COMMON.IOUNITS'
 569       external ufparm
 570       integer uiparm(1),k
 571       double precision urparm(1)
 572       dimension x(maxvar),g(maxvar)
 573 c
 574 c
 575 c
 576 cbad      icg=mod(nf,2)+1
 577       icg=1
 578 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 579       if (nf-nfl+1) 20,30,40
 580    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 581 cd      print *,20
 582       if (nf.eq.0) return
 583       goto 40
 584    30 continue
 585 cd      print *,30
 586       k=0
 587       do i=1,nres-1
 588         do j=1,3
 589           k=k+1
 590           dc(j,i)=x(k)
 591         enddo
 592       enddo
 593       do i=2,nres-1
 594         if (ialph(i,1).gt.0) then
 595         do j=1,3
 596           k=k+1
 597           dc(j,i+nres)=x(k)
 598         enddo
 599         endif
 600       enddo
 601       call chainbuild_cart
 602
 603 C
 604 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 605 C
 606    40 call cartgrad
 607 cd      print *,40
 608       k=0
 609       do i=1,nres-1
 610         do j=1,3
 611           k=k+1
 612           g(k)=gcart(j,i)
 613         enddo
 614       enddo
 615       do i=2,nres-1
 616         if (ialph(i,1).gt.0) then
 617         do j=1,3
 618           k=k+1
 619           g(k)=gxcart(j,i)
 620         enddo
 621         endif
 622       enddo
 623
 624       return
 625       end