source/unres/src-5hdiag-tmp/q_measure3.F

   1       double precision function qwolynes(seg1,seg2,flag,seg3,seg4)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.INTERACT'
   7       include 'COMMON.VAR'
   8       include 'COMMON.QRESTR'
   9       integer i,j,jl,k,l,il,kl,nl,np,seg1,seg2,seg3,seg4,secseg
  10       integer nsep /3/
  11       double precision dist,qm
  12       double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
  13       logical lprn /.false./
  14       logical flag
  15       qq = 0.0d0
  16       nl=0
  17       do i=0,nres
  18         do j=1,3
  19           dqwol(j,i)=0.0d0
  20           dxqwol(j,i)=0.0d0
  21         enddo
  22       enddo
  23       if (lprn) then
  24       write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4,
  25      & " flag",flag
  26       call flush(iout)
  27       endif
  28       if (flag) then
  29         do il=seg1+nsep,seg2
  30           do jl=seg1,il-nsep
  31             nl=nl+1
  32             if (itype(il).ne.10) then
  33               ilnres=il+nres
  34             else
  35               ilnres=il
  36             endif
  37             if (itype(jl).ne.10) then
  38               jlnres=jl+nres
  39             else
  40               jlnres=jl
  41             endif
  42             qqijCM = qcontrib(il,jl,ilnres,jlnres)
  43             qq = qq+qqijCM
  44             if (lprn) then
  45               write (iout,*) "qqijCM",qqijCM
  46               call flush(iout)
  47             endif
  48           enddo
  49         enddo
  50         if (lprn) then
  51           write (iout,*) "nl",nl," qq",qq
  52           call flush(iout)
  53         endif
  54       else
  55         do il=seg1,seg2
  56           if((seg3-il).lt.3) then
  57              secseg=il+3
  58           else
  59              secseg=seg3
  60           endif
  61           do jl=secseg,seg4
  62             nl=nl+1
  63             if (itype(il).ne.10) then
  64               ilnres=il+nres
  65             else
  66               ilnres=il
  67             endif
  68             if (itype(jl).ne.10) then
  69               jlnres=jl+nres
  70             else
  71               jlnres=jl
  72             endif
  73             qqijCM = qcontrib(il,jl,ilnres,jlnres)
  74             qq = qq+qqijCM
  75             if (lprn) then
  76               write (iout,*) "qqijCM",qqijCM
  77               call flush(iout)
  78             endif
  79           enddo
  80         enddo
  81       endif
  82       qq = qq/nl
  83       qwolynes=1.0d0-qq
  84       do i=0,nres
  85         do j=1,3
  86           dqwol(j,i)=dqwol(j,i)/nl
  87           dxqwol(j,i)=dxqwol(j,i)/nl
  88         enddo
  89       enddo
  90       return
  91       end
  92 #ifdef QDEBUG
  93 c-------------------------------------------------------------------
  94       subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
  95       implicit real*8 (a-h,o-z)
  96       include 'DIMENSIONS'
  97       include 'COMMON.IOUNITS'
  98       include 'COMMON.CHAIN'
  99       include 'COMMON.INTERACT'
 100       include 'COMMON.VAR'
 101       include 'COMMON.QRESTR'
 102       integer seg1,seg2,seg3,seg4
 103       logical flag
 104       double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
 105      & qwolxan(3,0:maxres),q1,q2
 106       double precision delta /1.0d-7/
 107       write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4
 108       write(iout,*) "dQ/dc backbone "
 109        do i=0,nres
 110         write(iout,'(i5,3e15.5)') i, (dqwol(j,i),j=1,3)
 111       enddo
 112       write(iout,*) "dQ/dX side chain "
 113       do i=1,nres
 114             write(iout,'(i5,3e15.5)') i,(dxqwol(j,i),j=1,3)
 115       enddo
 116       do i=1,nres
 117         do j=1,3
 118           cdummy(j,i)=c(j,i)
 119           c(j,i)=c(j,i)-delta
 120           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 121           c(j,i)=cdummy(j,i)+delta
 122           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 123           qwolan(j,i)=0.5d0*(q2-q1)/delta
 124           c(j,i)=cdummy(j,i)
 125 c          write (iout,*) "i",i," j",j," q1",q1," a2",q2
 126         enddo
 127       enddo
 128       do i=1,nres
 129         do j=1,3
 130           cdummy(j,i+nres)=c(j,i+nres)
 131           c(j,i+nres)=c(j,i+nres)-delta
 132           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 133           c(j,i+nres)=cdummy(j,i+nres)+delta
 134           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 135           qwolxan(j,i)=0.5d0*(q2-q1)/delta
 136           c(j,i+nres)=cdummy(j,i+nres)
 137         enddo
 138       enddo
 139       write(iout,*) "Numerical Q cartesian gradients backbone: "
 140       do i=0,nres
 141         write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
 142       enddo
 143       write(iout,*) "Numerical Q cartesian gradients side-chain: "
 144       do i=0,nres
 145         write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
 146       enddo
 147       return
 148       end
 149 #endif
 150 c------------------------------------------------------------------------
 151       subroutine EconstrQ
 152 c     MD with umbrella_sampling using Wolyne's distance measure as a constraint
 153       implicit real*8 (a-h,o-z)
 154       include 'DIMENSIONS'
 155       include 'COMMON.CONTROL'
 156       include 'COMMON.VAR'
 157       include 'COMMON.QRESTR'
 158 #ifndef LANG0
 159       include 'COMMON.LANGEVIN'
 160 #else
 161       include 'COMMON.LANGEVIN.lang0'
 162 #endif
 163       include 'COMMON.CHAIN'
 164       include 'COMMON.DERIV'
 165       include 'COMMON.GEO'
 166       include 'COMMON.LOCAL'
 167       include 'COMMON.INTERACT'
 168       include 'COMMON.IOUNITS'
 169       include 'COMMON.NAMES'
 170       include 'COMMON.TIME1'
 171       double precision uzap1,uzap2,hm1,hm2,hmnum
 172       double precision ucdelan,dUcartan(3,0:MAXRES)
 173      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
 174      &  duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
 175       integer kstart,kend,lstart,lend,idummy
 176       double precision delta /1.0d-7/
 177       do i=0,nres
 178          do j=1,3
 179             duconst(j,i)=0.0d0
 180             dudconst(j,i)=0.0d0
 181             duxconst(j,i)=0.0d0
 182             dudxconst(j,i)=0.0d0
 183          enddo
 184       enddo
 185       Uconst=0.0d0
 186       do i=1,nfrag
 187          qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
 188      &   ,idummy,idummy)
 189          Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
 190 c Calculating the derivatives of Constraint energy with respect to Q
 191          Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),qinfrag(i,iset))
 192 c Calculating the derivatives of Q with respect to cartesian coordinates
 193          do ii=0,nres
 194             do j=1,3
 195                duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
 196                dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
 197             enddo
 198          enddo
 199 c      write (iout,*) "Calling qwol_num"
 200 c      call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.,idummy,idummy)
 201       enddo
 202 c      stop
 203       do i=1,npair
 204          kstart=ifrag(1,ipair(1,i,iset),iset)
 205          kend=ifrag(2,ipair(1,i,iset),iset)
 206          lstart=ifrag(1,ipair(2,i,iset),iset)
 207          lend=ifrag(2,ipair(2,i,iset),iset)
 208          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
 209          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
 210 c  Calculating dU/dQ
 211          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
 212 c Calculating dQ/dXi
 213          do ii=0,nres
 214             do j=1,3
 215                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
 216                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
 217             enddo
 218          enddo
 219       enddo
 220 c      write(iout,*) "Uconst inside subroutine ", Uconst
 221 c Transforming the gradients from Cs to dCs for the backbone
 222       do i=0,nres
 223          do j=i+1,nres
 224            do k=1,3
 225              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
 226            enddo
 227          enddo
 228       enddo
 229 c  Transforming the gradients from Cs to dCs for the side chains
 230       do i=1,nres
 231          do j=1,3
 232            dudxconst(j,i)=duxconst(j,i)
 233          enddo
 234       enddo
 235 c      write(iout,*) "dU/dc backbone "
 236 c       do ii=0,nres
 237 c        write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
 238 c      enddo
 239 c      write(iout,*) "dU/dX side chain "
 240 c      do ii=1,nres
 241 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
 242 c      enddo
 243 c      write(iout,*) "dU/ddc backbone "
 244 c       do ii=0,nres
 245 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
 246 c      enddo
 247 c      write(iout,*) "dU/ddX side chain "
 248 c      do ii=1,nres
 249 c            write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
 250 c      enddo
 251 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
 252 c      call dEconstrQ_num
 253       return
 254       end
 255 #ifdef QDEBUG
 256 c-----------------------------------------------------------------------
 257       subroutine dEconstrQ_num
 258 c Calculating numerical dUconst/ddc and dUconst/ddx
 259       implicit real*8 (a-h,o-z)
 260       include 'DIMENSIONS'
 261       include 'COMMON.CONTROL'
 262       include 'COMMON.VAR'
 263       include 'COMMON.QRESTR'
 264 #ifndef LANG0
 265       include 'COMMON.LANGEVIN'
 266 #else
 267       include 'COMMON.LANGEVIN.lang0'
 268 #endif
 269       include 'COMMON.CHAIN'
 270       include 'COMMON.DERIV'
 271       include 'COMMON.GEO'
 272       include 'COMMON.LOCAL'
 273       include 'COMMON.INTERACT'
 274       include 'COMMON.IOUNITS'
 275       include 'COMMON.NAMES'
 276       include 'COMMON.TIME1'
 277       double precision uzap1,uzap2
 278       double precision dUcartan(3,0:MAXRES)
 279      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
 280       integer kstart,kend,lstart,lend,idummy
 281       double precision delta /1.0d-7/
 282 c     For the backbone
 283       do i=0,nres-1
 284          do j=1,3
 285             dUcartan(j,i)=0.0d0
 286             cdummy(j,i)=dc(j,i)
 287             dc(j,i)=dc(j,i)+delta
 288             call chainbuild_cart
 289             uzap2=0.0d0
 290             do ii=1,nfrag
 291                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 292      &           .true.,idummy,idummy)
 293                uzap2=uzap2+wfrag(ii,iset)*
 294      &                harmonic(qfrag(ii),qinfrag(ii,iset))
 295             enddo
 296             do ii=1,npair
 297                kstart=ifrag(1,ipair(1,ii,iset),iset)
 298                kend=ifrag(2,ipair(1,ii,iset),iset)
 299                lstart=ifrag(1,ipair(2,ii,iset),iset)
 300                lend=ifrag(2,ipair(2,ii,iset),iset)
 301                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 302                uzap2=uzap2+wpair(ii,iset)*
 303      &                harmonic(qpair(ii),qinpair(ii,iset))
 304             enddo
 305             dc(j,i)=cdummy(j,i)
 306             call chainbuild_cart
 307             uzap1=0.0d0
 308              do ii=1,nfrag
 309                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 310      &           .true.,idummy,idummy)
 311                uzap1=uzap1+wfrag(ii,iset)*
 312      &                harmonic(qfrag(ii),qinfrag(ii,iset))
 313             enddo
 314             do ii=1,npair
 315                kstart=ifrag(1,ipair(1,ii,iset),iset)
 316                kend=ifrag(2,ipair(1,ii,iset),iset)
 317                lstart=ifrag(1,ipair(2,ii,iset),iset)
 318                lend=ifrag(2,ipair(2,ii,iset),iset)
 319                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 320                uzap1=uzap1+wpair(ii,iset)*
 321      &                harmonic(qpair(ii),qinpair(ii,iset))
 322             enddo
 323             ducartan(j,i)=(uzap2-uzap1)/(delta)
 324          enddo
 325       enddo
 326 c Calculating numerical gradients for dU/ddx
 327       do i=0,nres-1
 328          do j=1,3
 329            duxcartan(j,i)=0.0d0
 330          enddo
 331          do j=1,3
 332             cdummy(j,i)=dc(j,i+nres)
 333             dc(j,i+nres)=dc(j,i+nres)+delta
 334             call chainbuild_cart
 335             uzap2=0.0d0
 336             do ii=1,nfrag
 337                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 338      &           .true.,idummy,idummy)
 339                uzap2=uzap2+wfrag(ii,iset)*
 340      &                harmonic(qfrag(ii),qinfrag(ii,iset))
 341             enddo
 342             do ii=1,npair
 343                kstart=ifrag(1,ipair(1,ii,iset),iset)
 344                kend=ifrag(2,ipair(1,ii,iset),iset)
 345                lstart=ifrag(1,ipair(2,ii,iset),iset)
 346                lend=ifrag(2,ipair(2,ii,iset),iset)
 347                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 348                uzap2=uzap2+wpair(ii,iset)*
 349      &                harmonic(qpair(ii),qinpair(ii,iset))
 350             enddo
 351             dc(j,i+nres)=cdummy(j,i)
 352             call chainbuild_cart
 353             uzap1=0.0d0
 354              do ii=1,nfrag
 355                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 356      &           .true.,idummy,idummy)
 357                uzap1=uzap1+wfrag(ii,iset)*
 358      &                 harmonic(qfrag(ii),qinfrag(ii,iset))
 359             enddo
 360             do ii=1,npair
 361                kstart=ifrag(1,ipair(1,ii,iset),iset)
 362                kend=ifrag(2,ipair(1,ii,iset),iset)
 363                lstart=ifrag(1,ipair(2,ii,iset),iset)
 364                lend=ifrag(2,ipair(2,ii,iset),iset)
 365                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 366                uzap1=uzap1+wpair(ii,iset)*
 367      &                harmonic(qpair(ii),qinpair(ii,iset))
 368             enddo
 369             duxcartan(j,i)=(uzap2-uzap1)/(delta)
 370          enddo
 371       enddo
 372       write(iout,*) "Numerical dUconst/ddc backbone "
 373       do ii=0,nres
 374         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
 375       enddo
 376       write(iout,*) "Numerical dUconst/ddx side-chain "
 377       do ii=1,nres
 378          write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
 379       enddo
 380       return
 381       end
 382 #endif
 383 c---------------------------------------------------------------------------
 384       double precision function qcontrib(il,jl,il1,jl1)
 385       implicit none
 386       include 'DIMENSIONS'
 387       include 'COMMON.IOUNITS'
 388       include 'COMMON.CHAIN'
 389       include 'COMMON.INTERACT'
 390       include 'COMMON.QRESTR'
 391       include 'COMMON.LOCAL'
 392       integer i,j,k,il,jl,il1,jl1,nd,itl,jtl
 393       double precision dist
 394       external dist
 395       double precision dij,dij1,d0ij,d0ij1,om1,om2,om12,om10,om20,om120
 396      &  ,fac,fac1,ddave,ssij,ddqij,d0ii1,d0jj1,rij,eom1,eom2,eom12
 397       double precision u(3),v(3),er(3),er0(3),dcosom1(3),dcosom2(3),
 398      &  aux1,aux2
 399       double precision scalar
 400       external scalar
 401       logical lprn /.false./
 402       if (lprn) write (iout,*) "il",il," jl",jl," il1",il1," jl1",jl1
 403       d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
 404      &           (cref(2,jl)-cref(2,il))**2+
 405      &           (cref(3,jl)-cref(3,il))**2)
 406       dij=dist(il,jl)
 407       dij1=dist(il1,jl1)
 408       do i=1,3
 409         er(i)=(c(i,jl1)-c(i,il1))/dij1
 410       enddo
 411       do i=1,3
 412         er0(i)=cref(i,jl1)-cref(i,il1)
 413       enddo
 414       d0ij1=dsqrt(scalar(er0,er0))
 415       do i=1,3
 416         er0(i)=er0(i)/d0ij1
 417       enddo
 418       if (il.ne.il1 .or. jl.ne.jl1) then
 419         ddave=0.5d0*((dij-d0ij)**2+(dij1-d0ij1)**2)
 420         nd=2
 421       else
 422         ddave=(dij-d0ij)**2
 423         nd=1
 424       endif
 425       if (il.ne.il1) then
 426         do i=1,3
 427           u(i)=cref(i,il1)-cref(i,il)
 428         enddo
 429         d0ii1=dsqrt(scalar(u,u))
 430         do i=1,3
 431           u(i)=u(i)/d0ii1
 432         enddo
 433         if (lprn) then
 434         write (iout,*) "u",(u(i),i=1,3)
 435         write (iout,*) "er0",(er0(i),i=1,3)
 436         om10=scalar(er0,u)
 437         om1=scalar(er,dc_norm(1,il1))
 438         write (iout,*) "om10",om10," om1",om1
 439         endif
 440       else
 441         om1=0.0d0
 442         om10=0.0d0
 443       endif
 444       if (jl.ne.jl1) then
 445         do i=1,3
 446           v(i)=cref(i,jl1)-cref(i,jl)
 447         enddo
 448         d0jj1=dsqrt(scalar(v,v))
 449         do i=1,3
 450           v(i)=v(i)/d0jj1
 451         enddo
 452         if (lprn) then
 453         write (iout,*) "v",(v(i),i=1,3)
 454         write (iout,*) "er0",(er0(i),i=1,3)
 455         om20=scalar(er,v)
 456         om2=scalar(er,dc_norm(1,jl1))
 457         write (iout,*) "om20",om20," om2",om2
 458         endif
 459       else
 460         om2=0.0d0
 461         om20=0.0d0
 462       endif
 463       if (il.ne.il1 .and. jl.ne.jl1) then
 464         om120=scalar(u,v)
 465         om12=scalar(dc_norm(1,il1),dc_norm(1,jl1))
 466       else
 467         om12=0.0d0
 468         om120=0.0d0
 469       endif
 470       if (lprn) then
 471         write (iout,*) "il",il," jl",jl,itype(il),itype(jl)
 472         write (iout,*)"d0ij",d0ij," om10",om10," om20",om20,
 473      &   " om120",om120,
 474      &  " dij",dij," om1",om1," om2",om2," om12",om12
 475         call flush(iout)
 476       endif
 477       ssij = 16.0d0/(d0ij*d0ij)
 478       qcontrib = dexp(-0.5d0*(ddave*ssij+((om1-om10)**2
 479      &                       +(om2-om20)**2+(om12-om120)**2)))
 480       if (lprn) write (iout,*) "ssij",ssij," qcontrib",qcontrib
 481 c      qcontrib = dexp(-0.5d0*(ddave*ssij)+(om1-om10)**2+(om2-om20)**2)
 482 c      qcontrib = dexp(-0.5d0*(ddave*ssij))
 483 c Compute gradient - radial component
 484       fac1 = qcontrib*ssij/nd
 485       fac = fac1*(dij-d0ij)/dij
 486       do k=1,3
 487         ddqij = (c(k,il)-c(k,jl))*fac
 488         dqwol(k,il)=dqwol(k,il)+ddqij
 489         dqwol(k,jl)=dqwol(k,jl)-ddqij
 490       enddo
 491       if (il1.ne.il .or. jl1.ne.jl) then
 492         fac = fac1*(dij1-d0ij1)/dij1
 493         do k=1,3
 494           ddqij = (c(k,il1)-c(k,jl1))*fac
 495           if (il1.ne.il) then
 496             dxqwol(k,il)=dxqwol(k,il)+ddqij
 497           else
 498             dqwol(k,il)=dqwol(k,il)+ddqij
 499           endif
 500           if (jl1.ne.jl) then
 501             dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 502           else
 503             dqwol(k,jl)=dqwol(k,jl)-ddqij
 504           endif
 505         enddo
 506       endif
 507 c      return
 508 c Orientational contributions
 509       rij=1.0d0/dij1
 510       eom1=qcontrib*(om1-om10)
 511       eom2=qcontrib*(om2-om20)
 512       eom12=qcontrib*(om12-om120)
 513       do k=1,3
 514         dcosom1(k)=rij*(dc_norm(k,il1)-om1*er(k))
 515         dcosom2(k)=rij*(dc_norm(k,jl1)-om2*er(k))
 516       enddo
 517       do k=1,3
 518         ddqij=eom1*dcosom1(k)+eom2*dcosom2(k)
 519         aux1=(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
 520      &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
 521         aux2=(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
 522      &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
 523         dqwol(k,il)=dqwol(k,il)-ddqij-aux1
 524         dqwol(k,jl)=dqwol(k,jl)+ddqij-aux2
 525         dxqwol(k,il)=dxqwol(k,il)-ddqij+aux1
 526 c     &            +(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
 527 c     &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
 528         dxqwol(k,jl)=dxqwol(k,jl)+ddqij+aux2
 529 c     &            +(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
 530 c     &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
 531       enddo
 532       return
 533       end