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