source/unres/src_CSA/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 cmd          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 cmd          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 cmd          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 cmd          call fricmat_mult(z,d_a_tmp)
 230         else if (iorder.eq.10) then
 231 cmd          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 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