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