source/analysis/src-M-prop/ave_property.F

   1       subroutine ave_property(ncon,*)
   2 ! construct the conformational ensembles at REMD temperatures
   3       implicit none
   4       include "DIMENSIONS"
   5       include "sizesclu.dat"
   6 #ifdef MPI
   7       include "mpif.h"
   8       include "COMMON.MPI"
   9       integer ierror,errcode,status(MPI_STATUS_SIZE)
  10 #endif
  11       include "COMMON.IOUNITS"
  12       include "COMMON.FREE"
  13       include "COMMON.FFIELD"
  14       include "COMMON.INTERACT"
  15       include "COMMON.SBRIDGE"
  16       include "COMMON.CHAIN"
  17       include "COMMON.CLUSTER"
  18       include "COMMON.PROPERTY"
  19       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  20      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
  21       double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc,
  22      &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  23      &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
  24      &      evdw_t,beta,w
  25       integer i,ii,ik,iproc,iscor,j,k,l,ib,ncon
  26       double precision sumprob,eini,efree,rmsdev
  27       character*80 bxname
  28       character*2 licz1
  29       character*5 ctemper
  30       integer ilen
  31       external ilen
  32       real*8 Fmin,Fmin_all
  33       double precision energia(0:max_ene)
  34       double precision qfree(600),ave_rms_chain(600),ave_assoc(600),
  35      & ave_hbond_intra(600),ave_hbond_inter(600),ave_assoc_hbond(600)
  36 #ifdef MPI
  37       double precision qfree_part(600),ave_rms_chain_part(600),
  38      & ave_assoc_part(600),ave_hbond_intra_part(600),
  39      & ave_hbond_inter_part(600),ave_assoc_hbond_part(600)
  40 c      write (iout,*) me," indstart",indstart(me)," indend",indend(me)
  41       call daread_ccoords(indstart(me),indend(me))
  42 #endif
  43
  44 #ifdef MPI
  45       do i=1,scount(me)
  46 #else
  47       do i=1,ncon
  48 #endif
  49         do j=1,nres
  50           do k=1,3
  51             c(k,j)=allcart(k,j,i)
  52             c(k,j+nres)=allcart(k,j+nres,i)
  53           enddo
  54         enddo
  55         do k=1,3
  56           c(k,nres+1)=c(k,1)
  57           c(k,nres+nres)=c(k,nres)
  58         enddo
  59         nss=nss_all(i)
  60         do j=1,nss
  61           ihpb(j)=ihpb_all(j,i)
  62           jhpb(j)=jhpb_all(j,i)
  63         enddo
  64         call int_from_cart1(.false.)
  65         do k=1,6
  66            ft(k)=1.0d0
  67         enddo
  68         call etotal(energia(0),fT)
  69 #ifdef DEBUG
  70         call enerprint(energia,ft)
  71 #endif
  72         do k=1,max_ene
  73           enetb(k,i)=energia(k)
  74         enddo
  75         call rmsd_intra(i)
  76         call monomer_contact(i)
  77         call hbonds_intra_inter(i)
  78 #ifdef DEBUG
  79         write (iout,*) "jcon",i
  80         write (iout,*) (rms_intrachain(j,i),j=1,nchain)
  81         write (iout,*) (npept_cont_interchain(j,i),j=1,nchain)
  82         write (iout,*) (npept_cont_intrachain(j,i),j=1,nchain)
  83         write (iout,*) (iassoc(j,i),j=1,nchain)
  84         call flush(iout)
  85 #endif
  86       enddo
  87 #ifdef DEBUG
  88       do i=1,ncon
  89         write (iout,*) "jcon",i
  90         write (iout,*) (rms_intrachain(j,i),j=1,nchain)
  91         write (iout,*) (npept_cont_interchain(j,i),j=1,nchain)
  92         write (iout,*) (npept_cont_intrachain(j,i),j=1,nchain)
  93         write (iout,*) (iassoc(j,i),j=1,nchain)
  94       enddo
  95 #endif
  96       do ib=200,500
  97
  98       beta=1.0d0/(1.987d-3*ib)
  99       Fmin=1.0d10
 100 #ifdef MPI
 101       do i=1,scount(me)
 102         ii=i+indstart(me)-1
 103 #else
 104       do i=1,ncon
 105         ii=i
 106 #endif
 107         evdw=enetb(1,i)
 108 #ifdef SCP14
 109         evdw2_14=enetb(17,i)
 110         evdw2=enetb(2,i)+evdw2_14
 111 #else
 112         evdw2=enetb(2,i)
 113         evdw2_14=0.0d0
 114 #endif
 115 #ifdef SPLITELE
 116         ees=enetb(3,i)
 117         evdw1=enetb(16,i)
 118 #else
 119         ees=enetb(3,i)
 120         evdw1=0.0d0
 121 #endif
 122         ecorr=enetb(4,i)
 123         ecorr5=enetb(5,i)
 124         ecorr6=enetb(6,i)
 125         eel_loc=enetb(7,i)
 126         eello_turn3=enetb(8,i)
 127         eello_turn4=enetb(9,i)
 128         eturn6=enetb(10,i)
 129         ebe=enetb(11,i)
 130         escloc=enetb(12,i)
 131         etors=enetb(13,i)
 132         etors_d=enetb(14,i)
 133         ehpb=enetb(15,i)
 134         estr=enetb(18,i)
 135         esccor=enetb(19,i)
 136         edihcnstr=enetb(20,i)
 137         evdw_t=enetb(21,i)
 138 #ifdef DEBUG
 139         write (iout,*) "evdw", wsc, evdw,evdw_t
 140         write (iout,*) "evdw2", wscp, evdw2
 141         write (iout,*) "welec", ft(1),welec,ees
 142         write (iout,*) "evdw1", wvdwpp,evdw1
 143         write (iout,*) "ebe" ebe,wang
 144 #endif
 145 #ifdef DEBUG
 146         write (iout,*) "fdim calc", i,ii,
 147      &   totfree(i),entfac(ii),Fdimless(i)
 148 #endif
 149         temper=ib+0.0d0
 150         if (rescale_mode.eq.1) then
 151           quot=1.0d0/(T0*beta*1.987D-3)
 152           quotl=1.0d0
 153           kfacl=1.0d0
 154           do l=1,5
 155             quotl1=quotl
 156             quotl=quotl*quot
 157             kfacl=kfacl*kfac
 158             fT(l)=kfacl/(kfacl-1.0d0+quotl)
 159           enddo
 160 #if defined(FUNCTH)
 161           ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 162      &            320.0d0
 163           ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 164           ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 165      &           /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 166 #elif defined(FUNCT)
 167           fT(6)=betaT/T0
 168           ftprim(6)=1.0d0/T0
 169           ftbis(6)=0.0d0
 170 #else
 171           fT(6)=1.0d0
 172           ftprim(6)=0.0d0
 173           ftbis(6)=0.0d0
 174 #endif
 175         else if (rescale_mode.eq.2) then
 176           quot=1.0d0/(T0*beta*1.987D-3)
 177           quotl=1.0d0
 178           do l=1,5
 179             quotl=quotl*quot
 180             fT(l)=1.12692801104297249644d0/
 181      &        dlog(dexp(quotl)+dexp(-quotl))
 182           enddo
 183 #if defined(FUNCTH)
 184           ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 185      &            320.0d0
 186           ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 187           ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 188      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 189 #elif defined(FUNCT)
 190           fT(6)=betaT/T0
 191           ftprim(6)=1.0d0/T0
 192           ftbis(6)=0.0d0
 193 #else
 194           fT(6)=1.0d0
 195           ftprim(6)=0.0d0
 196           ftbis(6)=0.0d0
 197 #endif
 198         endif
 199 c        write (iout,*) "ft",(ft(k),k=1,6)
 200 #ifdef SPLITELE
 201         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+welec*ft(1)*ees
 202      &   +wvdwpp*evdw1
 203      &   +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
 204      &   +wstrain*ehpb+nss*ebr+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
 205      &   +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
 206      &   +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
 207      &   +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
 208      &   +wbond*estr+wsccor*ft(1)*esccor
 209 #else
 210         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 211      &   +welec*ft(1)*(ees+evdw1)
 212      &   +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
 213      &   +wstrain*ehpb+nss*ebr+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
 214      &   +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
 215      &   +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
 216      &   +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
 217      &   +wbond*estr+wsccor*ft(1)*esccor
 218 #endif
 219 #ifdef DEBUG
 220       write (iout,10) (evdw+ft(6)*evdw_t),wsc,evdw2,wscp,ees,
 221      &  welec*ft(1),evdw1,
 222      &  wvdwpp,
 223      &  estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*ft(1),
 224      &  etors_d,wtor_d*ft(2),ehpb,wstrain,
 225      &  ecorr,wcorr*ft(3),ecorr5,wcorr5*ft(4),ecorr6,wcorr6*ft(5),
 226      &  eel_loc,wel_loc*ft(2),eello_turn3,wturn3*ft(2),
 227      &  eello_turn4,wturn4*ft(3),eturn6,wturn6*ft(5),
 228      &  esccor,wsccor*ft(1),edihcnstr,ebr*nss,etot
 229    10 format (/'Virtual-chain energies:'//
 230      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
 231      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
 232      & 'EES=   ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p elec)'/
 233      & 'EVDWPP=',1pE16.6,' WEIGHT=',1pD16.6,' (p-p VDW)'/
 234      & 'ESTR=  ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/
 235      & 'EBE=   ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/
 236      & 'ESC=   ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
 237      & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
 238      & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
 239      & 'EHBP=  ',1pE16.6,' WEIGHT=',1pD16.6,
 240      & ' (SS bridges & dist. cnstr.)'/
 241      & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
 242      & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
 243      & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
 244      & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/
 245      & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/
 246      & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/
 247      & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
 248      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
 249      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
 250      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
 251      & 'ETOT=  ',1pE16.6,' (total)')
 252 #endif
 253         totfree(i)=beta*etot+entfac(ii)
 254 #ifdef DEBUG
 255         write (iout,*) "ib",ib," beta",beta," i",i," ii",ii,
 256      &    " etot",etot,entfac(ii),totfree(i)
 257 #endif
 258         if (totfree(i).lt.Fmin) Fmin=totfree(i)
 259       enddo
 260 #ifdef MPI
 261       call MPI_Allreduce(Fmin,Fmin_all,1,
 262      &   MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,IERROR)
 263 #else
 264       Fmin_all=Fmin
 265 #endif
 266       qfree_part(ib)=0.0d0
 267       ave_rms_chain_part(ib)=0.0d0
 268       ave_assoc_part(ib)=0.0d0
 269       ave_hbond_intra_part(ib)=0.0d0
 270       ave_hbond_inter_part(ib)=0.0d0
 271       ave_assoc_hbond_part(ib)=0.0d0
 272 #ifdef MPI
 273       do i=1,scount(me)
 274         w=dexp(-totfree(i)+Fmin_all)
 275 c        write (iout,*) "i",i," w",w
 276         qfree_part(ib)=qfree_part(ib)+w
 277         w=w/nchain
 278         do j=1,nchain
 279           ave_rms_chain_part(ib)=ave_rms_chain_part(ib)
 280      &           +w*rms_intrachain(j,i)
 281           ave_assoc_part(ib)=ave_assoc_part(ib)+w*iassoc(j,i)
 282           ave_hbond_intra_part(ib)=ave_hbond_intra_part(ib)
 283      &           +w*npept_cont_intrachain(j,i)
 284           ave_hbond_inter_part(ib)=ave_hbond_inter_part(ib)
 285      &           +w*npept_cont_interchain(j,i)
 286           if (npept_cont_interchain(j,i).gt.0)
 287      &     ave_assoc_hbond_part(ib)=ave_assoc_hbond_part(ib)+w
 288         enddo
 289       enddo
 290
 291       enddo
 292
 293       call MPI_Reduce(qfree_part(1),qfree(1),600,
 294      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR)
 295       call MPI_Reduce(ave_rms_chain_part(1),ave_rms_chain(1),600,
 296      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR)
 297       call MPI_Reduce(ave_assoc_part(1),ave_assoc(1),600,
 298      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR)
 299       call MPI_Reduce(ave_hbond_intra_part(1),ave_hbond_intra(1),600,
 300      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR)
 301       call MPI_Reduce(ave_hbond_inter_part(1),ave_hbond_inter(1),600,
 302      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR)
 303       call MPI_Reduce(ave_assoc_hbond_part(1),ave_assoc_hbond(1),600,
 304      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR)
 305 #else
 306       qfree(ib)=0.0d0
 307       ave_rms_chain(ib)=0.0d0
 308       ave_assoc(ib)=0.0d0
 309       ave_hbond_intra(ib)=0.0d0
 310       ave_hbond_inter(ib)=0.0d0
 311       ave_assoc_hbond(ib)=0.0d0
 312       do i=1,ncon
 313         w=dexp(-totfree(i)+Fmin_all)
 314         qfree(ib)=qfree(ib)+w
 315         w=w/nchain
 316         do j=1,nchain
 317           ave_rms_chain(ib)=ave_rms_chain(ib)+w*rms_intrachain(j,i)
 318           ave_assoc(ib)=ave_assoc(ib)+w*iassoc(j,i)
 319           ave_hbond_intra(ib)=ave_hbond_intra(ib)
 320      &           +w*npept_cont_intrachain(j,i)
 321           ave_hbond_inter(ib)=ave_hbond_inter(ib)
 322      &           +w*npept_cont_interchain(j,i)
 323           if (npept_cont_interchain(j,i).gt.0)
 324      &     ave_assoc_hbond(ib)=ave_assoc_hbond(ib)+w
 325         enddo
 326       enddo
 327
 328       enddo
 329
 330 #endif
 331 #ifdef MPI
 332       if (me.eq.Master) then
 333 #endif
 334       do ib=200,500
 335         write (istat,'(i5,5f10.5)') ib,ave_rms_chain(ib)/qfree(ib),
 336      &   ave_assoc(ib)/qfree(ib),ave_assoc_hbond(ib)/qfree(ib),
 337      &   ave_hbond_intra(ib)/qfree(ib),ave_hbond_inter(ib)/qfree(ib)
 338       enddo
 339       call flush(istat)
 340 #ifdef MPI
 341       endif
 342 #endif
 343       return
 344       end
 345 !--------------------------------------------------
 346       subroutine mysort1(n, x, ipermut)
 347       implicit none
 348       integer i,j,imax,ipm,n
 349       real x(n)
 350       integer ipermut(n)
 351       real xtemp
 352       do i=1,n
 353         xtemp=x(i)
 354         imax=i
 355         do j=i+1,n
 356           if (x(j).lt.xtemp) then
 357             imax=j
 358             xtemp=x(j)
 359           endif
 360         enddo
 361         x(imax)=x(i)
 362         x(i)=xtemp
 363         ipm=ipermut(imax)
 364         ipermut(imax)=ipermut(i)
 365         ipermut(i)=ipm
 366       enddo
 367       return
 368       end
 369 !--------------------------------------------------
 370       subroutine rmsd_intra(jcon)
 371       implicit none
 372       include "DIMENSIONS"
 373       include "sizesclu.dat"
 374       include "COMMON.CONTROL"
 375       include "COMMON.IOUNITS"
 376       include "COMMON.FREE"
 377       include "COMMON.FFIELD"
 378       include "COMMON.INTERACT"
 379       include "COMMON.SBRIDGE"
 380       include "COMMON.CHAIN"
 381       include "COMMON.CLUSTER"
 382       include "COMMON.PROPERTY"
 383       integer i,ii,j,k
 384       double precision xx(3,maxres),yy(3,maxres)
 385       double precision rms
 386       integer jcon
 387       logical non_conv
 388       double precision przes(3),obrot(3,3)
 389       if (lside) then
 390         do i=1,nchain
 391         ii=0
 392         do k=ichain(1,i),ichain(2,i)
 393           ii=ii+1
 394           do j=1,3
 395             xx(j,ii)=allcart(j,i,jcon)
 396             yy(j,ii)=cref(j,i,1)
 397           enddo
 398         enddo
 399         do k=ichain(1,i),ichain(2,i)
 400 c          if (itype(i).ne.10) then
 401             ii=ii+1
 402             do j=1,3
 403               xx(j,ii)=allcart(j,i+nres,jcon)
 404               yy(j,ii)=cref(j,i+nres,1)
 405             enddo
 406 c          endif
 407         enddo
 408         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 409         rms_intrachain(i,jcon)=dsqrt(rms)
 410         enddo
 411       else
 412         do i=nnt,nct
 413           do j=1,3
 414             xx(j,i)=allcart(j,i,jcon)
 415           enddo
 416         enddo
 417 #ifdef DEBUG
 418         write (iout,*) "Conformation",jcon
 419 #endif
 420         do i=1,nchain
 421 #ifdef DEBUG
 422           write (iout,*) "chain",i
 423 #endif
 424           do j=1,ichain(2,i)-ichain(1,i)+1
 425 #ifdef DEBUG
 426             write (iout,'(i5,3f10.5,5x,3f10.5)') j,
 427      &        (xx(k,ichain(1,i)+j-1),k=1,3),
 428      &        (cref(k,ichain(1,i)+j-1,1),k=1,3)
 429 #endif
 430           enddo
 431           call fitsq(rms,xx(1,ichain(1,i)),cref(1,ichain(1,i),1),
 432      &        ichain(2,i)-ichain(1,i)+1,przes,obrot,non_conv)
 433           rms_intrachain(i,jcon)=dsqrt(rms)
 434 #ifdef DEBUG
 435           write (iout,*) "rms",rms
 436 #endif
 437         enddo
 438       endif
 439       return
 440       end
 441 !--------------------------------------------------
 442       subroutine hbonds_intra_inter(jcon)
 443       implicit none
 444       include "DIMENSIONS"
 445       include "sizesclu.dat"
 446       include "COMMON.IOUNITS"
 447       include "COMMON.FREE"
 448       include "COMMON.FFIELD"
 449       include "COMMON.INTERACT"
 450       include "COMMON.SBRIDGE"
 451       include "COMMON.CHAIN"
 452       include "COMMON.CLUSTER"
 453       include "COMMON.PROPERTY"
 454       integer ncont,icont(2,maxcont)
 455       integer i,j,jcon,nr1,nr2
 456       do i=1,nchain
 457         npept_cont_interchain(i,jcon)=0
 458         npept_cont_intrachain(i,jcon)=0
 459       enddo
 460       do i=1,nres
 461         do j=1,3
 462           c(j,i)=allcart(j,i,jcon)
 463         enddo
 464       enddo
 465       call elecont(.false.,ncont,icont,nnt,nct)
 466       do i=1,ncont
 467         nr1=nres_chain(icont(1,i))
 468         nr2=nres_chain(icont(2,i))
 469         if (nr1.eq.nr2) then
 470           npept_cont_intrachain(nr1,jcon)=
 471      &     npept_cont_intrachain(nr1,jcon)+1
 472         else
 473           npept_cont_interchain(nr1,jcon)=
 474      &     npept_cont_interchain(nr1,jcon)+1
 475           npept_cont_interchain(nr2,jcon)=
 476      &     npept_cont_interchain(nr2,jcon)+1
 477         endif
 478       enddo
 479       return
 480       end
 481 !--------------------------------------------------
 482       subroutine monomer_contact(jcon)
 483       implicit none
 484       include "DIMENSIONS"
 485       include "sizesclu.dat"
 486       include "COMMON.IOUNITS"
 487       include "COMMON.FREE"
 488       include "COMMON.FFIELD"
 489       include "COMMON.INTERACT"
 490       include "COMMON.SBRIDGE"
 491       include "COMMON.CHAIN"
 492       include "COMMON.CLUSTER"
 493       include "COMMON.PROPERTY"
 494       integer ncont,icont(2,maxcont)
 495       integer i,j,ic1,ic2,jcon
 496       do i=1,nchain
 497         iassoc(i,jcon)=0
 498       enddo
 499       do i=1,nres
 500         do j=1,3
 501           c(j,i)=allcart(j,i,jcon)
 502         enddo
 503       enddo
 504       call contact(.false.,ncont,icont,nnt,nct)
 505       do i=1,ncont
 506         ic1=nres_chain(icont(1,i))
 507         ic2=nres_chain(icont(2,i))
 508         if (ic1.ne.ic2) then
 509           iassoc(ic1,jcon)=1
 510           iassoc(ic2,jcon)=1
 511         endif
 512       enddo
 513       return
 514       end