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