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