source/unres/src_MD-M-SAXS-homology/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 c      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,*) "in func x"
 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,*) "upon exit from func"
 274 c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 275 c      write (iout,*) 'f=',f
 276 c       jjj=0
 277 c     endif
 278       return
 279       end
 280 ************************************************************************
 281       subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
 282       implicit real*8 (a-h,o-z)
 283       include 'DIMENSIONS'
 284       include 'COMMON.DERIV'
 285       include 'COMMON.IOUNITS'
 286       include 'COMMON.GEO'
 287       common /chuju/ jjj
 288       double precision energia(0:n_ene)
 289       integer jjj
 290       double precision ufparm
 291       external ufparm
 292       integer uiparm(1)
 293       real*8 urparm(1)
 294       dimension x(maxvar)
 295 c     if (jjj.gt.0) then
 296 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 297 c     endif
 298       nfl=nf
 299       icg=mod(nf,2)+1
 300       call var_to_geom_restr(n,x)
 301       call zerograd
 302       call chainbuild_extconf
 303 cd    write (iout,*) 'ETOTAL called from FUNC'
 304       call etotal(energia(0))
 305       call sum_gradient
 306       f=energia(0)
 307 c     if (jjj.gt.0) then
 308 c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
 309 c       write (iout,*) 'f=',etot
 310 c       jjj=0
 311 c     endif
 312       return
 313       end
 314 c-------------------------------------------------------
 315       subroutine x2xx(x,xx,n)
 316       implicit real*8 (a-h,o-z)
 317       include 'DIMENSIONS'
 318       include 'COMMON.VAR'
 319       include 'COMMON.CHAIN'
 320       include 'COMMON.INTERACT'
 321       double precision xx(maxvar),x(maxvar)
 322
 323       do i=1,nvar
 324         varall(i)=x(i)
 325       enddo
 326
 327       ig=0
 328       igall=0
 329       do i=4,nres
 330         igall=igall+1
 331         if (mask_phi(i).eq.1) then
 332           ig=ig+1
 333           xx(ig)=x(igall)
 334         endif
 335       enddo
 336
 337       do i=3,nres
 338         igall=igall+1
 339         if (mask_theta(i).eq.1) then
 340           ig=ig+1
 341           xx(ig)=x(igall)
 342         endif
 343       enddo
 344
 345       do ij=1,2
 346       do i=2,nres-1
 347         if (itype(i).ne.10) then
 348           igall=igall+1
 349           if (mask_side(i).eq.1) then
 350             ig=ig+1
 351             xx(ig)=x(igall)
 352           endif
 353         endif
 354       enddo
 355       enddo
 356
 357       n=ig
 358
 359       return
 360       end
 361
 362       subroutine xx2x(x,xx)
 363       implicit real*8 (a-h,o-z)
 364       include 'DIMENSIONS'
 365       include 'COMMON.VAR'
 366       include 'COMMON.CHAIN'
 367       include 'COMMON.INTERACT'
 368       double precision xx(maxvar),x(maxvar)
 369
 370       do i=1,nvar
 371         x(i)=varall(i)
 372       enddo
 373
 374       ig=0
 375       igall=0
 376       do i=4,nres
 377         igall=igall+1
 378         if (mask_phi(i).eq.1) then
 379           ig=ig+1
 380           x(igall)=xx(ig)
 381         endif
 382       enddo
 383
 384       do i=3,nres
 385         igall=igall+1
 386         if (mask_theta(i).eq.1) then
 387           ig=ig+1
 388           x(igall)=xx(ig)
 389         endif
 390       enddo
 391
 392       do ij=1,2
 393       do i=2,nres-1
 394         if (itype(i).ne.10) then
 395           igall=igall+1
 396           if (mask_side(i).eq.1) then
 397             ig=ig+1
 398             x(igall)=xx(ig)
 399           endif
 400         endif
 401       enddo
 402       enddo
 403
 404       return
 405       end
 406
 407 c----------------------------------------------------------
 408       subroutine minim_dc(etot,iretcode,nfun)
 409       implicit real*8 (a-h,o-z)
 410       include 'DIMENSIONS'
 411       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 412 #ifdef MPI
 413       include 'mpif.h'
 414 #endif
 415       include 'COMMON.CONTROL'
 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       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 c-----
 491 c      write (iout,*) "checkgrad before SUMSL"
 492 c      icheckgrad=1
 493 c      aincr=1.0d-7
 494 c      call exec_checkgrad
 495 c-----
 496
 497       call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
 498 c-----
 499 c      write (iout,*) "checkgrad after SUMSL"
 500 c      call exec_checkgrad
 501 c-----
 502
 503       k=0
 504       do i=1,nres-1
 505         do j=1,3
 506           k=k+1
 507           dc(j,i)=x(k)
 508         enddo
 509       enddo
 510       do i=2,nres-1
 511         if (ialph(i,1).gt.0) then
 512         do j=1,3
 513           k=k+1
 514           dc(j,i+nres)=x(k)
 515         enddo
 516         endif
 517       enddo
 518       call chainbuild_cart
 519
 520 cd      call zerograd
 521 cd      nf=0
 522 cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
 523 cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
 524 cd
 525 cd      do i=1,k
 526 cd       x(i)=x(i)+1.0D-5
 527 cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 528 cd       x(i)=x(i)-1.0D-5
 529 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 530 cd      enddo
 531
 532       etot=v(10)
 533       iretcode=iv(1)
 534       nfun=iv(6)
 535       return
 536       end
 537
 538       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 539       implicit real*8 (a-h,o-z)
 540       include 'DIMENSIONS'
 541 #ifdef MPI
 542       include 'mpif.h'
 543 #endif
 544       include 'COMMON.SETUP'
 545       include 'COMMON.DERIV'
 546       include 'COMMON.IOUNITS'
 547       include 'COMMON.GEO'
 548       include 'COMMON.CHAIN'
 549       include 'COMMON.VAR'
 550       double precision energia(0:n_ene)
 551       double precision ufparm
 552       external ufparm
 553       integer uiparm(1)
 554       real*8 urparm(1)
 555       dimension x(maxvar)
 556       nfl=nf
 557 cbad      icg=mod(nf,2)+1
 558       icg=1
 559
 560       k=0
 561       do i=1,nres-1
 562         do j=1,3
 563           k=k+1
 564           dc(j,i)=x(k)
 565         enddo
 566       enddo
 567       do i=2,nres-1
 568         if (ialph(i,1).gt.0) then
 569         do j=1,3
 570           k=k+1
 571           dc(j,i+nres)=x(k)
 572         enddo
 573         endif
 574       enddo
 575       call chainbuild_cart
 576
 577       call zerograd
 578       call etotal(energia(0))
 579       f=energia(0)
 580
 581 cd      print *,'func_dc ',nf,nfl,f
 582
 583       return
 584       end
 585
 586       subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
 587       implicit real*8 (a-h,o-z)
 588       include 'DIMENSIONS'
 589 #ifdef MPI
 590       include 'mpif.h'
 591 #endif
 592       include 'COMMON.SETUP'
 593       include 'COMMON.CHAIN'
 594       include 'COMMON.DERIV'
 595       include 'COMMON.VAR'
 596       include 'COMMON.INTERACT'
 597       include 'COMMON.FFIELD'
 598       include 'COMMON.MD'
 599       include 'COMMON.IOUNITS'
 600       external ufparm
 601       integer uiparm(1),k
 602       double precision urparm(1)
 603       dimension x(maxvar),g(maxvar)
 604 c
 605 c
 606 c
 607 cbad      icg=mod(nf,2)+1
 608       icg=1
 609 cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
 610       if (nf-nfl+1) 20,30,40
 611    20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
 612 cd      print *,20
 613       if (nf.eq.0) return
 614       goto 40
 615    30 continue
 616 cd      print *,30
 617       k=0
 618       do i=1,nres-1
 619         do j=1,3
 620           k=k+1
 621           dc(j,i)=x(k)
 622         enddo
 623       enddo
 624       do i=2,nres-1
 625         if (ialph(i,1).gt.0) then
 626         do j=1,3
 627           k=k+1
 628           dc(j,i+nres)=x(k)
 629         enddo
 630         endif
 631       enddo
 632       call chainbuild_cart
 633
 634 C
 635 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 636 C
 637    40 call cartgrad
 638 cd      print *,40
 639       k=0
 640       do i=1,nres-1
 641         do j=1,3
 642           k=k+1
 643           g(k)=gcart(j,i)
 644         enddo
 645       enddo
 646       do i=2,nres-1
 647         if (ialph(i,1).gt.0) then
 648         do j=1,3
 649           k=k+1
 650           g(k)=gxcart(j,i)
 651         enddo
 652         endif
 653       enddo
 654
 655       return
 656       end