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