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       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       IF (mask_r) THEN
  88        call x2xx(x,xx,nvar_restr)
  89        call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
  90      &                    iv,liv,lv,v,idum,rdum,fdum)
  91        call xx2x(x,xx)
  92       ELSE
  93        call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
  94       ENDIF
  95       etot=v(10)
  96       iretcode=iv(1)
  97 cd    print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
  98 cd    write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
  99 c     call intout
 100 c     change=reduce(x)
 101       call var_to_geom(nvar,x)
 102 c     if (change) then
 103 c       write (iout,'(a)') 'Reduction worked, minimizing again...'
 104 c     else
 105 c       not_done=.false.
 106 c     endif
 107       call chainbuild
 108 c     call etotal(energia(0))
 109 c     etot=energia(0)
 110 c     call enerprint(energia(0))
 111       nfun=iv(6)
 112
 113 c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
 114
 115 c     ENDDO ! NOT_DONE
 116
 117       return
 118       end
 119 #ifdef MPI
 120 c----------------------------------------------------------------------------
 121       subroutine ergastulum
 122       implicit real*8 (a-h,o-z)
 123       include 'DIMENSIONS'
 124 #ifdef MPI
 125       include "mpif.h"
 126 #endif
 127       include 'COMMON.SETUP'
 128       include 'COMMON.DERIV'
 129       include 'COMMON.VAR'
 130       include 'COMMON.IOUNITS'
 131       include 'COMMON.FFIELD'
 132       include 'COMMON.INTERACT'
 133       include 'COMMON.MD'
 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         time00=MPI_Wtime()
 146         call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
 147         time_Bcast=time_Bcast+MPI_Wtime()-time00
 148         if (icall.gt.4 .and. iorder.ge.0)
 149      &   time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
 150        icall=icall+1
 151 c      write (*,*)
 152 c     & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
 153         if (iorder.eq.0) then
 154           call zerograd
 155           call etotal(edum)
 156 c          write (2,*) "After etotal"
 157 c          write (2,*) "dimen",dimen," dimen3",dimen3
 158 c          call flush(2)
 159         else if (iorder.eq.2) then
 160           call zerograd
 161           call etotal_short(edum)
 162 c          write (2,*) "After etotal_short"
 163 c          write (2,*) "dimen",dimen," dimen3",dimen3
 164 c          call flush(2)
 165         else if (iorder.eq.3) then
 166           call zerograd
 167           call etotal_long(edum)
 168 c          write (2,*) "After etotal_long"
 169 c          write (2,*) "dimen",dimen," dimen3",dimen3
 170 c          call flush(2)
 171         else if (iorder.eq.1) then
 172           call sum_gradient
 173 c          write (2,*) "After sum_gradient"
 174 c          write (2,*) "dimen",dimen," dimen3",dimen3
 175 c          call flush(2)
 176         else if (iorder.eq.4) then
 177           call ginv_mult(z,d_a_tmp)
 178         else if (iorder.eq.5) then
 179 c Setup MD things for a slave
 180           dimen=(nct-nnt+1)+nside
 181           dimen1=(nct-nnt)+(nct-nnt+1)
 182           dimen3=dimen*3
 183 c          write (2,*) "dimen",dimen," dimen3",dimen3
 184 c          call flush(2)
 185           call int_bounds(dimen,igmult_start,igmult_end)
 186           igmult_start=igmult_start-1
 187           call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 188      &     ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 189            my_ng_count=igmult_end-igmult_start
 190           call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 191      &     MPI_INTEGER,FG_COMM,IERROR)
 192 c          write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 193 c          write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 194           myginv_ng_count=maxres2*my_ng_count
 195 c          write (2,*) "igmult_start",igmult_start," igmult_end",
 196 c     &     igmult_end," my_ng_count",my_ng_count
 197 c          call flush(2)
 198           call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 199      &     nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 200           call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 201      &     nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 202 c          write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 203 c          write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 204 c          call flush(2)
 205 c          call MPI_Barrier(FG_COMM,IERROR)
 206           time00=MPI_Wtime()
 207           call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 208      &     nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 209      &     myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 210 #ifdef TIMING
 211           time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 212 #endif
 213           do i=1,dimen
 214             do j=1,2*my_ng_count
 215               ginv(j,i)=gcopy(i,j)
 216             enddo
 217           enddo
 218 c          write (2,*) "dimen",dimen," dimen3",dimen3
 219 c          write (2,*) "End MD setup"
 220 c          call flush(2)
 221 c           write (iout,*) "My chunk of ginv_block"
 222 c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
 223         else if (iorder.eq.6) then
 224           call int_from_cart1(.false.)
 225         else if (iorder.eq.7) then
 226           call chainbuild_cart
 227         else if (iorder.eq.8) then
 228           call intcartderiv
 229         else if (iorder.eq.9) then
 230           call fricmat_mult(z,d_a_tmp)
 231         else if (iorder.eq.10) then
 232           call setup_fricmat
 233         endif
 234       enddo
 235       write (*,*) 'Processor',fg_rank,' CG group',kolor,
 236      &  ' absolute rank',myrank,' leves ERGASTULUM.'
 237       write(*,*)'Processor',fg_rank,' wait times for respective orders',
 238      & (' order[',i,']',time_order(i),i=0,10)
 239       return
 240       end
 241 #endif
 242 ************************************************************************
 243       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
 244       implicit real*8 (a-h,o-z)
 245       include 'DIMENSIONS'
 246       include 'COMMON.DERIV'
 247       include 'COMMON.IOUNITS'
 248       include 'COMMON.GEO'
 249       common /chuju/ jjj
 250       double precision energia(0:n_ene)
 251       integer jjj
 252       double precision ufparm
 253       external ufparm
 254       integer uiparm(1)
 255       real*8 urparm(1)
 256       dimension x(maxvar)
 257 c     if (jjj.gt.0) then
 258 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 259 c     endif
 260       nfl=nf
 261       icg=mod(nf,2)+1
 262 cd      print *,'func',nf,nfl,icg
 263       call var_to_geom(n,x)
 264       call zerograd
 265       call chainbuild
 266 cd    write (iout,*) 'ETOTAL called from FUNC'
 267       call etotal(energia(0))
 268       call sum_gradient
 269       f=energia(0)
 270 c     if (jjj.gt.0) then
 271 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 272 c       write (iout,*) 'f=',etot
 273 c       jjj=0
 274 c     endif
 275       return
 276       end
 277 ************************************************************************
 278       subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
 279       implicit real*8 (a-h,o-z)
 280       include 'DIMENSIONS'
 281       include 'COMMON.DERIV'
 282       include 'COMMON.IOUNITS'
 283       include 'COMMON.GEO'
 284       common /chuju/ jjj
 285       double precision energia(0:n_ene)
 286       integer jjj
 287       double precision ufparm
 288       external ufparm
 289       integer uiparm(1)
 290       real*8 urparm(1)
 291       dimension x(maxvar)
 292 c     if (jjj.gt.0) then
 293 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 294 c     endif
 295       nfl=nf
 296       icg=mod(nf,2)+1
 297       call var_to_geom_restr(n,x)
 298       call zerograd
 299       call chainbuild
 300 cd    write (iout,*) 'ETOTAL called from FUNC'
 301       call etotal(energia(0))
 302       call sum_gradient
 303       f=energia(0)
 304 c     if (jjj.gt.0) then
 305 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 306 c       write (iout,*) 'f=',etot
 307 c       jjj=0
 308 c     endif
 309       return
 310       end
 311 c-------------------------------------------------------
 312       subroutine x2xx(x,xx,n)
 313       implicit real*8 (a-h,o-z)
 314       include 'DIMENSIONS'
 315       include 'COMMON.VAR'
 316       include 'COMMON.CHAIN'
 317       include 'COMMON.INTERACT'
 318       double precision xx(maxvar),x(maxvar)
 319
 320       do i=1,nvar
 321         varall(i)=x(i)
 322       enddo
 323
 324       ig=0
 325       igall=0
 326       do i=4,nres
 327         igall=igall+1
 328         if (mask_phi(i).eq.1) then
 329           ig=ig+1
 330           xx(ig)=x(igall)
 331         endif
 332       enddo
 333
 334       do i=3,nres
 335         igall=igall+1
 336         if (mask_theta(i).eq.1) then
 337           ig=ig+1
 338           xx(ig)=x(igall)
 339         endif
 340       enddo
 341
 342       do ij=1,2
 343       do i=2,nres-1
 344         if (itype(i).ne.10) then
 345           igall=igall+1
 346           if (mask_side(i).eq.1) then
 347             ig=ig+1
 348             xx(ig)=x(igall)
 349           endif
 350         endif
 351       enddo
 352       enddo
 353
 354       n=ig
 355
 356       return
 357       end
 358
 359       subroutine xx2x(x,xx)
 360       implicit real*8 (a-h,o-z)
 361       include 'DIMENSIONS'
 362       include 'COMMON.VAR'
 363       include 'COMMON.CHAIN'
 364       include 'COMMON.INTERACT'
 365       double precision xx(maxvar),x(maxvar)
 366
 367       do i=1,nvar
 368         x(i)=varall(i)
 369       enddo
 370
 371       ig=0
 372       igall=0
 373       do i=4,nres
 374         igall=igall+1
 375         if (mask_phi(i).eq.1) then
 376           ig=ig+1
 377           x(igall)=xx(ig)
 378         endif
 379       enddo
 380
 381       do i=3,nres
 382         igall=igall+1
 383         if (mask_theta(i).eq.1) then
 384           ig=ig+1
 385           x(igall)=xx(ig)
 386         endif
 387       enddo
 388
 389       do ij=1,2
 390       do i=2,nres-1
 391         if (itype(i).ne.10) then
 392           igall=igall+1
 393           if (mask_side(i).eq.1) then
 394             ig=ig+1
 395             x(igall)=xx(ig)
 396           endif
 397         endif
 398       enddo
 399       enddo
 400
 401       return
 402       end
 403
 404 c----------------------------------------------------------
 405       subroutine minim_dc(etot,iretcode,nfun)
 406       implicit real*8 (a-h,o-z)
 407       include 'DIMENSIONS'
 408       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 409 #ifdef MPI
 410       include 'mpif.h'
 411 #endif
 412       include 'COMMON.SETUP'
 413       include 'COMMON.IOUNITS'
 414       include 'COMMON.VAR'
 415       include 'COMMON.GEO'
 416       include 'COMMON.MINIM'
 417       include 'COMMON.CHAIN'
 418       dimension iv(liv)
 419       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
 420 c      common /przechowalnia/ v
 421
 422       double precision energia(0:n_ene)
 423       external func_dc,grad_dc,fdum
 424       logical not_done,change,reduce
 425       double precision g(maxvar),f1
 426
 427       call deflt(2,iv,liv,lv,v)
 428 * 12 means fresh start, dont call deflt
 429       iv(1)=12
 430 * max num of fun calls
 431       if (maxfun.eq.0) maxfun=500
 432       iv(17)=maxfun
 433 * max num of iterations
 434       if (maxmin.eq.0) maxmin=1000
 435       iv(18)=maxmin
 436 * controls output
 437       iv(19)=2
 438 * selects output unit
 439       iv(21)=0
 440       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 441 * 1 means to print out result
 442       iv(22)=print_min_res
 443 * 1 means to print out summary stats
 444       iv(23)=print_min_stat
 445 * 1 means to print initial x and d
 446       iv(24)=print_min_ini
 447 * min val for v(radfac) default is 0.1
 448       v(24)=0.1D0
 449 * max val for v(radfac) default is 4.0
 450       v(25)=2.0D0
 451 c     v(25)=4.0D0
 452 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 453 * the sumsl default is 0.1
 454       v(26)=0.1D0
 455 * false conv if (act fnctn decrease) .lt. v(34)
 456 * the sumsl default is 100*machep
 457       v(34)=v(34)/100.0D0
 458 * absolute convergence
 459       if (tolf.eq.0.0D0) tolf=1.0D-4
 460       v(31)=tolf
 461 * relative convergence
 462       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 463       v(32)=rtolf
 464 * controls initial step size
 465        v(35)=1.0D-1
 466 * large vals of d correspond to small components of step
 467       do i=1,6*nres
 468         d(i)=1.0D-1
 469       enddo
 470
 471       k=0
 472       do i=1,nres-1
 473         do j=1,3
 474           k=k+1
 475           x(k)=dc(j,i)
 476         enddo
 477       enddo
 478       do i=2,nres-1
 479         if (ialph(i,1).gt.0) then
 480         do j=1,3
 481           k=k+1
 482           x(k)=dc(j,i+nres)
 483         enddo
 484         endif
 485       enddo
 486
 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