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