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_extconf
 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_extconf
 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.CONTROL'
 414       include 'COMMON.SETUP'
 415       include 'COMMON.IOUNITS'
 416       include 'COMMON.VAR'
 417       include 'COMMON.GEO'
 418       include 'COMMON.MINIM'
 419       include 'COMMON.CHAIN'
 420       dimension iv(liv)
 421       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
 422       common /przechowalnia/ v
 423
 424       double precision energia(0:n_ene)
 425       external func_dc,grad_dc,fdum
 426       logical not_done,change,reduce
 427       double precision g(maxvar),f1
 428
 429       call deflt(2,iv,liv,lv,v)
 430 * 12 means fresh start, dont call deflt
 431       iv(1)=12
 432 * max num of fun calls
 433       if (maxfun.eq.0) maxfun=500
 434       iv(17)=maxfun
 435 * max num of iterations
 436       if (maxmin.eq.0) maxmin=1000
 437       iv(18)=maxmin
 438 * controls output
 439       iv(19)=2
 440 * selects output unit
 441       iv(21)=0
 442       if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
 443 * 1 means to print out result
 444       iv(22)=print_min_res
 445 * 1 means to print out summary stats
 446       iv(23)=print_min_stat
 447 * 1 means to print initial x and d
 448       iv(24)=print_min_ini
 449 * min val for v(radfac) default is 0.1
 450       v(24)=0.1D0
 451 * max val for v(radfac) default is 4.0
 452       v(25)=2.0D0
 453 c     v(25)=4.0D0
 454 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 455 * the sumsl default is 0.1
 456       v(26)=0.1D0
 457 * false conv if (act fnctn decrease) .lt. v(34)
 458 * the sumsl default is 100*machep
 459       v(34)=v(34)/100.0D0
 460 * absolute convergence
 461       if (tolf.eq.0.0D0) tolf=1.0D-4
 462       v(31)=tolf
 463 * relative convergence
 464       if (rtolf.eq.0.0D0) rtolf=1.0D-4
 465       v(32)=rtolf
 466 * controls initial step size
 467        v(35)=1.0D-1
 468 * large vals of d correspond to small components of step
 469       do i=1,6*nres
 470         d(i)=1.0D-1
 471       enddo
 472
 473       k=0
 474       do i=1,nres-1
 475         do j=1,3
 476           k=k+1
 477           x(k)=dc(j,i)
 478         enddo
 479       enddo
 480       do i=2,nres-1
 481         if (ialph(i,1).gt.0) then
 482         do j=1,3
 483           k=k+1
 484           x(k)=dc(j,i+nres)
 485         enddo
 486         endif
 487       enddo
 488 c-----
 489 c      write (iout,*) "checkgrad before SUMSL"
 490 c      icheckgrad=1
 491 c      aincr=1.0d-7
 492 c      call exec_checkgrad
 493 c-----
 494
 495       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 496 c-----
 497 c      write (iout,*) "checkgrad after SUMSL"
 498 c      call exec_checkgrad
 499 c-----
 500
 501       k=0
 502       do i=1,nres-1
 503         do j=1,3
 504           k=k+1
 505           dc(j,i)=x(k)
 506         enddo
 507       enddo
 508       do i=2,nres-1
 509         if (ialph(i,1).gt.0) then
 510         do j=1,3
 511           k=k+1
 512           dc(j,i+nres)=x(k)
 513         enddo
 514         endif
 515       enddo
 516       call chainbuild_cart
 517
 518 cd      call zerograd
 519 cd      nf=0
 520 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 521 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 522 cd
 523 cd      do i=1,k
 524 cd       x(i)=x(i)+1.0D-5
 525 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 526 cd       x(i)=x(i)-1.0D-5
 527 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 528 cd      enddo
 529
 530       etot=v(10)
 531       iretcode=iv(1)
 532       nfun=iv(6)
 533       return
 534       end
 535
 536       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 537       implicit real*8 (a-h,o-z)
 538       include 'DIMENSIONS'
 539 #ifdef MPI
 540       include 'mpif.h'
 541 #endif
 542       include 'COMMON.SETUP'
 543       include 'COMMON.DERIV'
 544       include 'COMMON.IOUNITS'
 545       include 'COMMON.GEO'
 546       include 'COMMON.CHAIN'
 547       include 'COMMON.VAR'
 548       double precision energia(0:n_ene)
 549       double precision ufparm
 550       external ufparm
 551       integer uiparm(1)
 552       real*8 urparm(1)
 553       dimension x(maxvar)
 554       nfl=nf
 555 cbad      icg=mod(nf,2)+1
 556       icg=1
 557
 558       k=0
 559       do i=1,nres-1
 560         do j=1,3
 561           k=k+1
 562           dc(j,i)=x(k)
 563         enddo
 564       enddo
 565       do i=2,nres-1
 566         if (ialph(i,1).gt.0) then
 567         do j=1,3
 568           k=k+1
 569           dc(j,i+nres)=x(k)
 570         enddo
 571         endif
 572       enddo
 573       call chainbuild_cart
 574
 575       call zerograd
 576       call etotal(energia(0))
 577       f=energia(0)
 578
 579 cd      print *,'func_dc ',nf,nfl,f
 580
 581       return
 582       end
 583
 584       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 585       implicit real*8 (a-h,o-z)
 586       include 'DIMENSIONS'
 587 #ifdef MPI
 588       include 'mpif.h'
 589 #endif
 590       include 'COMMON.SETUP'
 591       include 'COMMON.CHAIN'
 592       include 'COMMON.DERIV'
 593       include 'COMMON.VAR'
 594       include 'COMMON.INTERACT'
 595       include 'COMMON.FFIELD'
 596       include 'COMMON.MD'
 597       include 'COMMON.IOUNITS'
 598       external ufparm
 599       integer uiparm(1),k
 600       double precision urparm(1)
 601       dimension x(maxvar),g(maxvar)
 602 c
 603 c
 604 c
 605 cbad      icg=mod(nf,2)+1
 606       icg=1
 607 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 608       if (nf-nfl+1) 20,30,40
 609    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 610 cd      print *,20
 611       if (nf.eq.0) return
 612       goto 40
 613    30 continue
 614 cd      print *,30
 615       k=0
 616       do i=1,nres-1
 617         do j=1,3
 618           k=k+1
 619           dc(j,i)=x(k)
 620         enddo
 621       enddo
 622       do i=2,nres-1
 623         if (ialph(i,1).gt.0) then
 624         do j=1,3
 625           k=k+1
 626           dc(j,i+nres)=x(k)
 627         enddo
 628         endif
 629       enddo
 630       call chainbuild_cart
 631
 632 C
 633 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 634 C
 635    40 call cartgrad
 636 cd      print *,40
 637       k=0
 638       do i=1,nres-1
 639         do j=1,3
 640           k=k+1
 641           g(k)=gcart(j,i)
 642         enddo
 643       enddo
 644       do i=2,nres-1
 645         if (ialph(i,1).gt.0) then
 646         do j=1,3
 647           k=k+1
 648           g(k)=gxcart(j,i)
 649         enddo
 650         endif
 651       enddo
 652
 653       return
 654       end