source/unres/src_MD-M/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       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       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 sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 487
 488       k=0
 489       do i=1,nres-1
 490         do j=1,3
 491           k=k+1
 492           dc(j,i)=x(k)
 493         enddo
 494       enddo
 495       do i=2,nres-1
 496         if (ialph(i,1).gt.0) then
 497         do j=1,3
 498           k=k+1
 499           dc(j,i+nres)=x(k)
 500         enddo
 501         endif
 502       enddo
 503       call chainbuild_cart
 504
 505 cd      call zerograd
 506 cd      nf=0
 507 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 508 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 509 cd
 510 cd      do i=1,k
 511 cd       x(i)=x(i)+1.0D-5
 512 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 513 cd       x(i)=x(i)-1.0D-5
 514 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 515 cd      enddo
 516
 517       etot=v(10)
 518       iretcode=iv(1)
 519       nfun=iv(6)
 520       return
 521       end
 522
 523       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 524       implicit real*8 (a-h,o-z)
 525       include 'DIMENSIONS'
 526 #ifdef MPI
 527       include 'mpif.h'
 528 #endif
 529       include 'COMMON.SETUP'
 530       include 'COMMON.DERIV'
 531       include 'COMMON.IOUNITS'
 532       include 'COMMON.GEO'
 533       include 'COMMON.CHAIN'
 534       include 'COMMON.VAR'
 535       double precision energia(0:n_ene)
 536       double precision ufparm
 537       external ufparm
 538       integer uiparm(1)
 539       real*8 urparm(1)
 540       dimension x(maxvar)
 541       nfl=nf
 542 cbad      icg=mod(nf,2)+1
 543       icg=1
 544
 545       k=0
 546       do i=1,nres-1
 547         do j=1,3
 548           k=k+1
 549           dc(j,i)=x(k)
 550         enddo
 551       enddo
 552       do i=2,nres-1
 553         if (ialph(i,1).gt.0) then
 554         do j=1,3
 555           k=k+1
 556           dc(j,i+nres)=x(k)
 557         enddo
 558         endif
 559       enddo
 560       call chainbuild_cart
 561
 562       call zerograd
 563       call etotal(energia(0))
 564       f=energia(0)
 565
 566 cd      print *,'func_dc ',nf,nfl,f
 567
 568       return
 569       end
 570
 571       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 572       implicit real*8 (a-h,o-z)
 573       include 'DIMENSIONS'
 574 #ifdef MPI
 575       include 'mpif.h'
 576 #endif
 577       include 'COMMON.SETUP'
 578       include 'COMMON.CHAIN'
 579       include 'COMMON.DERIV'
 580       include 'COMMON.VAR'
 581       include 'COMMON.INTERACT'
 582       include 'COMMON.FFIELD'
 583       include 'COMMON.MD'
 584       include 'COMMON.IOUNITS'
 585       external ufparm
 586       integer uiparm(1),k
 587       double precision urparm(1)
 588       dimension x(maxvar),g(maxvar)
 589 c
 590 c
 591 c
 592 cbad      icg=mod(nf,2)+1
 593       icg=1
 594 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 595       if (nf-nfl+1) 20,30,40
 596    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 597 cd      print *,20
 598       if (nf.eq.0) return
 599       goto 40
 600    30 continue
 601 cd      print *,30
 602       k=0
 603       do i=1,nres-1
 604         do j=1,3
 605           k=k+1
 606           dc(j,i)=x(k)
 607         enddo
 608       enddo
 609       do i=2,nres-1
 610         if (ialph(i,1).gt.0) then
 611         do j=1,3
 612           k=k+1
 613           dc(j,i+nres)=x(k)
 614         enddo
 615         endif
 616       enddo
 617       call chainbuild_cart
 618
 619 C
 620 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 621 C
 622    40 call cartgrad
 623 cd      print *,40
 624       k=0
 625       do i=1,nres-1
 626         do j=1,3
 627           k=k+1
 628           g(k)=gcart(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           g(k)=gxcart(j,i)
 636         enddo
 637         endif
 638       enddo
 639
 640       return
 641       end