subroutine read_constr_homology include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.HOMRESTR' include 'COMMON.HOMOLOGY' c c For new homol impl c include 'COMMON.VAR' c include 'include_unres/COMMON.VAR' c c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d, c & dist_cut c common /przechowalnia/ odl_temp(maxres,maxres,max_template), c & sigma_odl_temp(maxres,maxres,max_template) character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp integer idomain(max_template,maxres) logical lprn /.true./ integer ilen external ilen logical liiflag c c FP - Nov. 2014 Temporary specifications for new vars c double precision rescore_tmp,x12,y12,z12,rescore2_tmp, & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 double precision, dimension (max_template,maxres) :: rescore3 character*24 tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints c not using pre-computed ones stored in files model_ki_{dist,angle} c FP (Nov., 2014) c ----------------------------------------------------------------- c c c Alternative: reading from input call card_concat(controlcard,.true.) call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0) call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0) call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) call readi(controlcard,"IHSET",ihset,1) if (homol_nset.gt.1)then call card_concat(controlcard,.true.) read(controlcard,*) (waga_homology(i),i=1,homol_nset) c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then c write(iout,*) "iset homology_weight " c do i=1,homol_nset c write(iout,*) i,waga_homology(i) c enddo c endif iset=mod(kolor,homol_nset)+1 else iset=1 waga_homology(1)=1.0 endif c write(iout,*) "waga_homology(",iset,")",waga_homology(iset) cd write (iout,*) "nnt",nnt," nct",nct cd call flush(iout) lim_odl=0 lim_dih=0 c c New c lim_theta=0 lim_xx=0 c c Reading HM global scores (prob not required) c do i = nnt,nct do k=1,constr_homology idomain(k,i)=0 enddo enddo c open (4,file="HMscore") c do k=1,constr_homology c read (4,*,end=521) hmscore_tmp c hmscore(k)=hmscore_tmp ! Another transformation can be used c write(*,*) "Model", k, ":", hmscore(k) c enddo c521 continue ii=0 do i = nnt,nct-2 do j=i+2,nct ii=ii+1 ii_in_use(ii)=0 enddo enddo c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d if (read_homol_frag) then call read_klapaucjusz else do k=1,constr_homology read(inp,'(a)') pdbfile c Next stament causes error upon compilation (?) c if(me.eq.king.or. .not. out1file) c write (iout,'(2a)') 'PDB data will be read from file ', c & pdbfile(:ilen(pdbfile)) write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file', & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file', & pdbfile(:ilen(pdbfile)) stop 34 continue c print *,'Begin reading pdb data' c c Files containing res sim or local scores (former containing sigmas) c write(kic2,'(bz,i2.2)') k tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. if (read2sigma) then call readpdb_template(k) close(ipdbin) else call readpdb close(ipdbin) endif c call readpdb do i=1,2*nres do j=1,3 crefjlee(j,i)=c(j,i) enddo enddo #ifdef DEBUG do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo write (iout,*) "read_constr_homology: after reading pdb file" call flush(iout) #endif c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) do i=1,2*nres do j=1,3 c(j,i)=cref(j,i) c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo c c From read_dist_constr (commented out 25/11/2014 <-> res sim) c c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore open (ientin,file=tpl_k_rescore,status='old') if (nnt.gt.1) rescore(k,1)=0.0d0 do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, & rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp rescore3(k,i_tmp)=rescore3_tmp write(iout,'(a7,i5,3f10.5,i5)') "rescore", & i_tmp,rescore2_tmp,rescore_tmp, & rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) endif enddo 1401 continue close (ientin) if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 do j=i+2,nct x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", c & constr_homology ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 l_homo(k,ii)=.false. endif enddo enddo lim_odl=ii endif c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_dih,status='old') c do irec=1,maxres-3 ! loop for reading sigma_dih c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for? c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right? c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_dih(k,i+nnt-1) c enddo c1402 continue c close (ientin) do i = nnt+3,nct if (idomain(k,i).eq.0) then sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) ! right? c read (ientin,*) sigma_dih(k,i) ! original variant c write (iout,*) "dih(",k,i,") =",dih(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2), c & "rescore(",k,i-3,") =",rescore(k,i-3) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 c write (iout,*) "Raw sigmas for dihedral angle restraints" c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2)*rescore(k,i-3) ! right expression ? c Instead of res sim other local measure of b/b str reliability possible if (sigma_dih(k,i).ne.0) & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) enddo lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_theta,status='old') c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for? c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_theta(k,i+nnt-1) c enddo c1403 continue c close (ientin) do i = nnt+2,nct ! right? without parallel. c do i = i=1,nres ! alternative for bounds acc to readpdb? c do i=ithet_start,ithet_end ! with FG parallel. if (idomain(k,i).eq.0) then sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2) c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 if (sigma_theta(k,i).ne.0) & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for? c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_d(k,i+nnt-1) c enddo c1404 continue do i = nnt,nct ! right? without parallel. c do i=2,nres-1 ! alternative for bounds acc to readpdb? c do i=loc_start,loc_end ! with FG parallel. if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i) c sigma_d(k,i)=rescore(k,i) ! right expression ? sigma_d(k,i)=rescore3(k,i) ! right expression ? if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) c read (ientin,*) sigma_d(k,i) ! 1st variant enddo endif enddo c c remove distance restraints not used in any model from the list c shift data in all arrays c if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 if (ii_in_use(ii).eq.0.and.liiflag) then liiflag=.false. iistart=ii endif if (ii_in_use(ii).ne.0.and..not.liiflag.or. & .not.liiflag.and.ii.eq.lim_odl) then if (ii.eq.lim_odl) then iishift=ii-iistart+1 else iishift=ii-iistart endif liiflag=.true. do ki=iistart,lim_odl-iishift ires_homo(ki)=ires_homo(ki+iishift) jres_homo(ki)=jres_homo(ki+iishift) ii_in_use(ki)=ii_in_use(ki+iishift) do k=1,constr_homology odl(k,ki)=odl(k,ki+iishift) sigma_odl(k,ki)=sigma_odl(k,ki+iishift) l_homo(k,ki)=l_homo(k,ki+iishift) enddo enddo ii=ii-iishift lim_odl=lim_odl-iishift endif enddo enddo endif endif ! .not. klapaucjusz if (constr_homology.gt.0) call homology_partition if (constr_homology.gt.0) call init_int_table cd write (iout,*) "homology_partition: lim_theta= ",lim_theta, cd & "lim_xx=",lim_xx c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end c c Print restraints c if (.not.lprn) return cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" do ii=1,lim_odl write(iout,'(3i5,100(2f8.2,1x,l1,4x))') & ii,ires_homo(ii),jres_homo(ii), & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii), & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" do i=nnt+3,nct write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*dih(ki,i), & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" do i=nnt+2,nct write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*thetatpl(ki,i), & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) enddo c endif c ----------------------------------------------------------------- return end c---------------------------------------------------------------------- subroutine read_klapaucjusz include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.HOMOLOGY' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.HOMRESTR' character*256 fragfile integer ninclust(maxclust),inclust(max_template,maxclust), & nresclust(maxclust),iresclust(maxres,maxclust) character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp integer idomain(max_template,maxres) logical lprn /.true./ integer ilen external ilen logical liiflag c c double precision rescore_tmp,x12,y12,z12,rescore2_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 character*24 tpl_k_rescore c c For new homol impl c include 'COMMON.VAR' c call getenv("FRAGFILE",fragfile) open(ientin,file=fragfile,status="old",err=10) read(ientin,*) constr_homology,nclust l_homo = .false. sigma_theta=0.0 sigma_d=0.0 sigma_dih=0.0 c Read pdb files do k=1,constr_homology read(ientin,'(a)') pdbfile write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file', & pdbfile(:ilen(pdbfile)) stop 34 continue unres_pdb=.false. call readpdb_template(k) c do i=1,2*nres c do j=1,3 c chomo(j,i,k)=c(j,i) c enddo c enddo do i=1,nres rescore(k,i)=0.2d0 rescore2(k,i)=1.0d0 enddo enddo c Read clusters do i=1,nclust read(ientin,*) ninclust(i),nresclust(i) read(ientin,*) (inclust(k,i),k=1,ninclust(i)) read(ientin,*) (iresclust(k,i),k=1,nresclust(i)) enddo c c Loop over clusters c do l=1,nclust do ll = 1,ninclust(l) k = inclust(ll,l) do i=1,nres idomain(k,i)=0 enddo do i=1,nresclust(l) if (nnt.gt.1) then idomain(k,iresclust(i,l)+1) = 1 else idomain(k,iresclust(i,l)) = 1 endif enddo c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) c write (iout,*) "k",k do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo c call cartprint c write (iout,*) "read_klapaucjusz: calling int_from_cart" call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) c write (iout,*) "en" do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 do j=i+2,nct x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", c & constr_homology ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 c l_homo(k,ii)=.false. endif enddo enddo lim_odl=ii endif c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then do i = nnt+3,nct if (idomain(k,i).eq.0) then c sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), c & " sigma_dihed",sigma_dih(k,i) if (sigma_dih(k,i).ne.0) & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) enddo lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then do i = nnt+2,nct if (idomain(k,i).eq.0) then c sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 if (sigma_theta(k,i).ne.0) & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) enddo endif if (waga_d.gt.0.0d0) then do i = nnt,nct if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then c sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) sigma_d(k,i)=rescore(k,i) if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 enddo endif enddo ! l enddo ! ll c c remove distance restraints not used in any model from the list c shift data in all arrays c if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 if (ii_in_use(ii).eq.0.and.liiflag) then liiflag=.false. iistart=ii endif if (ii_in_use(ii).ne.0.and..not.liiflag.or. & .not.liiflag.and.ii.eq.lim_odl) then if (ii.eq.lim_odl) then iishift=ii-iistart+1 else iishift=ii-iistart endif liiflag=.true. do ki=iistart,lim_odl-iishift ires_homo(ki)=ires_homo(ki+iishift) jres_homo(ki)=jres_homo(ki+iishift) ii_in_use(ki)=ii_in_use(ki+iishift) do k=1,constr_homology odl(k,ki)=odl(k,ki+iishift) sigma_odl(k,ki)=sigma_odl(k,ki+iishift) l_homo(k,ki)=l_homo(k,ki+iishift) enddo enddo ii=ii-iishift lim_odl=lim_odl-iishift endif enddo enddo endif #ifdef DEBUG write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl do i=1,lim_odl write (iout,*) i,ires_homo(i),jres_homo(i) enddo #endif return 10 stop "Error in fragment file" end