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