2D replica exchange with homology constraints

diff --git a/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe
index 0929cef..b628330 100755 (executable)
Binary files a/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe and b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe differ

diff --git a/source/unres/src_MD/COMMON.CONTROL b/source/unres/src_MD/COMMON.CONTROL
index 9fce3c5..0202bae 100644 (file)
--- a/source/unres/src_MD/COMMON.CONTROL
+++ b/source/unres/src_MD/COMMON.CONTROL
@@ -1,5 +1,5 @@
       integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
-     & inprint,i2ndstr,mucadyn,constr_dist,constr_homology
+     & inprint,i2ndstr,mucadyn,constr_dist,constr_homology,homol_nset
       real*8 waga_dist, waga_angle
       logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
      &                 sideadd,lsecondary,read_cart,unres_pdb,
@@ -10,6 +10,7 @@
      & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
      & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
      & constr_dist,gnorm_check,gradout,split_ene,constr_homology,
-     & waga_dist, waga_angle
+     & homol_nset
+      common /homol/ waga_dist(maxprocs/20), waga_angle(maxprocs/20)
 C... minim = .true. means DO minimization.
 C... energy_dec = .true. means print energy decomposition matrix

diff --git a/source/unres/src_MD/MREMD.F b/source/unres/src_MD/MREMD.F
index be6af9c..9301122 100644 (file)
--- a/source/unres/src_MD/MREMD.F
+++ b/source/unres/src_MD/MREMD.F
@@ -70,6 +70,16 @@ cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
          enddo
       endif
 
+      if(homol_nset.gt.1) then
+         i_econstr=24
+         nset=homol_nset
+         do i=1,nset
+          mset(i)=1
+         enddo
+      endif
+
+      if(usampl) i_econstr=20
+
       k=0
       rep2i(k)=-1
       do il=1,max0(nset,1)
@@ -187,7 +197,7 @@ cd           write (*,*) me," After broadcast: file_exist",file_exist
               read (irest2,*) ndowna(0,il),
      &                    (ndowna(i,il),i=1,ndowna(0,il))
              enddo
-             if(usampl.or.hremd.gt.0) then
+             if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
               read (irest2,*)
               read (irest2,*) nset
               read (irest2,*) 
@@ -278,7 +288,7 @@ cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
          if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
 
        endif
-       if(usampl.or.hremd.gt.0) then
+       if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
           iset=i2set(me)
           if (hremd.gt.0) call set_hweights(iset)
           if(me.eq.king.or..not.out1file) 
@@ -572,7 +582,7 @@ c Variable time step algorithm.
               write (irest1,*) ndowna(0,il),
      &                   (ndowna(i,il),i=1,ndowna(0,il))
              enddo
-             if(usampl.or.hremd.gt.0) then
+             if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
               write (irest1,*) "nset"
               write (irest1,*) nset
               write (irest1,*) "mset"
@@ -599,7 +609,7 @@ c Variable time step algorithm.
            do i=1,2*nres
             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
            enddo
-           if(usampl.or.hremd.gt.0) then
+           if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
              write (irest2,*) iset
            endif
           close(irest2)
@@ -717,20 +727,31 @@ c     &             remd_t_bath,1,mpi_double_precision,king,
 c     &             CG_COMM,ierr)
            potEcomp(n_ene+1)=t_bath
            t_bath_old=t_bath
-           if (usampl) then
+           if (usampl.or.homol_nset.gt.1) then
              potEcomp(n_ene+2)=iset
+c             write(iout,*) potEcomp(n_ene+1),potEcomp(n_ene+2),iset,nset
              if (iset.lt.nset) then
                i_set_temp=iset
                iset=iset+1
-               call EconstrQ
-               potEcomp(n_ene+3)=Uconst
+               if (homol_nset.gt.1) then
+                call e_modeller(potEcomp(n_ene+3))
+c                write(iout,*) "iset+1",potEcomp(n_ene+3)
+               else
+                call EconstrQ
+                potEcomp(n_ene+3)=Uconst
+               endif
                iset=i_set_temp
              endif
              if (iset.gt.1) then
                i_set_temp=iset
                iset=iset-1
-               call EconstrQ
-               potEcomp(n_ene+4)=Uconst 
+               if (homol_nset.gt.1) then
+                call e_modeller(potEcomp(n_ene+4))
+c                write(iout,*) "iset-1",potEcomp(n_ene+4)
+               else
+                call EconstrQ
+                potEcomp(n_ene+4)=Uconst
+               endif
                iset=i_set_temp
              endif
            endif
@@ -778,9 +799,15 @@ ctime            call flush(iout)
 
 
           if (me.eq.king) then
+            write(iout,*) 
+     &     'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
             do i=1,nodes
                remd_t_bath(i)=remd_ene(n_ene+1,i)
                iremd_iset(i)=remd_ene(n_ene+2,i)
+               write(iout,'(i4,f10.3,f6.0,i3,2f10.3)') 
+     &                i,remd_ene(i_econstr,i),
+     &                remd_ene(n_ene+1,i),iremd_iset(i),
+     &                remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
             enddo
 #ifdef DEBUG
             if(lmuca) then
@@ -798,7 +825,7 @@ co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
             endif
 #endif
 c-------------------------------------           
-           IF(.not.usampl.and.hremd.eq.0) THEN
+           IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
 #ifdef DEBUG
             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
      &        " nodes",nodes
@@ -982,9 +1009,9 @@ c                call flush(iout)
            enddo
 cd           write (iout,*) "exchange completed"
 cd           call flush(iout) 
-        ELSEIF (usampl) THEN
+        ELSEIF (usampl.or.homol_nset.gt.1) THEN
           do ii=1,nodes  
-cd            write(iout,*) "########",ii
+c            write(iout,*) "########",ii
 
             i_temp=iran_num(1,nrep)
             i_mult=iran_num(1,remd_m(i_temp))
@@ -992,14 +1019,14 @@ cd            write(iout,*) "########",ii
             i_mset=iran_num(1,mset(i_iset))
             i=i_index(i_temp,i_mult,i_iset,i_mset)
 
-cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
+c            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
 
             if(t_exchange_only)then
              i_dir=1
             else
              i_dir=iran_num(1,3)
             endif
-cd            write(iout,*) "i_dir=",i_dir
+c            write(iout,*) "i_dir=",i_dir
 
             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then            
                
@@ -1016,10 +1043,11 @@ cd            write(iout,*) "i_dir=",i_dir
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
-               remd_ene(20,i)=remd_ene(n_ene+3,i)
-               remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+                
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
+               remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+               remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
 
             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
 
@@ -1028,16 +1056,16 @@ cd            write(iout,*) "i_dir=",i_dir
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
-               remd_ene(20,i)=remd_ene(n_ene+3,i)
-               remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
+               remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+               remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
 
             else
                goto 444 
             endif
  
-cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
+c            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
 ctime            call flush(iout)
 
 c Swap temperatures between conformations i and iex with recalculating the free energies
@@ -1050,33 +1078,39 @@ co     &          remd_t_bath(i)
               
               call sum_energy(remd_ene(0,iex),.false.)
               ene_iex_i=remd_ene(0,iex)
-cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
+cdebug
+c ERROR only makes sense for dir =1
+c              write (iout,*) "ene_iex_i",remd_ene(0,iex)
 c              call sum_energy(remd_ene(0,i),.false.)
-cd              write (iout,*) "ene_i_i",remd_ene(0,i)
+c              write (iout,*) "ene_i_i",remd_ene(0,i)
 c              write (iout,*) "rescaling weights with temperature",
 c     &          remd_t_bath(iex)
 c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
-c                write (iout,*) "ERROR: inconsistent energies:",i,
+c                write (iout,*) "ERROR: inconsistent energies i:",i,
 c     &            ene_i_i,remd_ene(0,i)
 c              endif
+cdebug_end
               call rescale_weights(remd_t_bath(iex))
               call sum_energy(remd_ene(0,i),.false.)
 cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
               ene_i_iex=remd_ene(0,i)
+cdebug
+c ERROR only makes sense for dir =1
 c              call sum_energy(remd_ene(0,iex),.false.)
 c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
-c                write (iout,*) "ERROR: inconsistent energies:",iex,
+c                write (iout,*) "ERROR: inconsistent energies iex:",iex,
 c     &            ene_iex_iex,remd_ene(0,iex)
 c              endif
-cd              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
+c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 c              write (iout,*) "i",i," iex",iex
-cd              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-cd     &           " ene_i_iex",ene_i_iex,
-cd     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
+c              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
+c     &           " ene_i_iex",ene_i_iex,
+c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
+cdebug_end
               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
               delta=-delta
-cd              write(iout,*) 'delta',delta
+c              write(iout,*) 'delta',delta
 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 c     &              (remd_t_bath(i)*remd_t_bath(iex))
@@ -1091,7 +1125,7 @@ c     &              (remd_t_bath(i)*remd_t_bath(iex))
      &          iremd_tot_usa(int(i2set(i-1)))=
      &                 iremd_tot_usa(int(i2set(i-1)))+1
               xxx=ran_number(0.0d0,1.0d0)
-cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
+c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
               if (delta .gt. xxx) then
                 tmp=remd_t_bath(i)       
                 remd_t_bath(i)=remd_t_bath(iex)
@@ -1127,8 +1161,8 @@ cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
               else
                remd_ene(0,iex)=ene_iex_iex
                remd_ene(0,i)=ene_i_i
-               remd_ene(20,iex)=econstr_temp_iex
-               remd_ene(20,i)=econstr_temp_i
+               remd_ene(i_econstr,iex)=econstr_temp_iex
+               remd_ene(i_econstr,i)=econstr_temp_i
               endif
 
 cd      do il=1,nset
@@ -1301,7 +1335,7 @@ c-------------------------------------
      &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
              enddo
 
-             if(usampl) then
+             if(usampl.or.homol_nset.gt.1) then
               do i=1,nset
                if(iremd_tot_usa(i).ne.0)
      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
@@ -1338,7 +1372,7 @@ cd         call flush(iout)
      &           CG_COMM,ierr) 
 cd         write (iout,*) "After scatter"
 cd         call flush(iout)
-         if(usampl.or.hremd.gt.0)
+         if(usampl.or.hremd.gt.0.or.homol_nset.gt.1)
      &    call mpi_scatter(iremd_iset,1,mpi_integer,
      &           iset,1,mpi_integer,king,
      &           CG_COMM,ierr) 
@@ -1511,7 +1545,7 @@ c-----------------------------------------------------------------------
            enddo
          enddo
 
-         if(usampl) then
+         if(usampl.or.homol_nset.gt.1) then
            call xdrfint_(ixdrf, nset, iret)
            do i=1,nset
              call xdrfint_(ixdrf,mset(i), iret)
@@ -1572,7 +1606,7 @@ c-----------------------------------------------------------------------
          enddo
 
 
-             if(usampl) then
+             if(usampl.or.homol_nset.gt.1) then
               call xdrfint(ixdrf, nset, iret)
               do i=1,nset
                 call xdrfint(ixdrf,mset(i), iret)
@@ -1954,7 +1988,7 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
          enddo
        
 
-           if(usampl) then
+           if(usampl.or.homol_nset.gt.1) then
 #ifdef AIX
              if(me.eq.king)then
               call xdrfint_(ixdrf, nset, iret)

diff --git a/source/unres/src_MD/cinfo.f b/source/unres/src_MD/cinfo.f
index b42001a..66c74d5 100644 (file)
--- a/source/unres/src_MD/cinfo.f
+++ b/source/unres/src_MD/cinfo.f
@@ -1,32 +1,29 @@
 C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-C 3 2 169
+C 3 2 191
       subroutine cinfo
       include 'COMMON.IOUNITS'
       write(iout,*)'++++ Compile info ++++'
-      write(iout,*)'Version 3.2 build 169'
-      write(iout,*)'compiled Mon Jul 21 16:29:49 2014'
-      write(iout,*)'compiled by jal47@matrix.chem.cornell.edu'
+      write(iout,*)'Version 3.2 build 191'
+      write(iout,*)'compiled Fri Nov  7 14:25:21 2014'
+      write(iout,*)'compiled by czarek@piasek4'
       write(iout,*)'OS name:    Linux '
-      write(iout,*)'OS release: 2.6.34.9-69.fc13.x86_64 '
+      write(iout,*)'OS release: 3.2.0-70-generic '
       write(iout,*)'OS version:',
-     & ' #1 SMP Tue May 3 09:23:03 UTC 2011 '
+     & ' #105-Ubuntu SMP Wed Sep 24 19:49:16 UTC 2014 '
       write(iout,*)'flags:'
       write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...'
-      write(iout,*)'FC = ifort'
-      write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include'
-      write(iout,*)'FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)...'
-      write(iout,*)'FFLAGS2 = -c -w -g -O0 -I$(INSTALL_DIR)/include'
-      write(iout,*)'FFLAGS3 = -c -w -O3 -mp'
-      write(iout,*)'FFLAGSE = -c -w -O3 -ipo -ipo_obj  -opt_report ...'
-      write(iout,*)'CC = cc'
-      write(iout,*)'CFLAGS = -DLINUX -DPGI -c'
-      write(iout,*)'OPT =  -O3 -ip -w'
+      write(iout,*)'FC= ifort'
+      write(iout,*)'OPT =  -O3 -ip '
+      write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include '
+      write(iout,*)'FFLAGS1 = -c  -g -CA -CB -I$(INSTALL_DIR)/inclu...'
+      write(iout,*)'FFLAGS2 = -c  -g -O0 -I$(INSTALL_DIR)/include  '
+      write(iout,*)'FFLAGSE = -c  -O3 -ipo  -opt_report -I$(INSTALL...'
       write(iout,*)'LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdr...'
       write(iout,*)'ARCH = LINUX'
       write(iout,*)'PP = /lib/cpp -P'
       write(iout,*)'object = unres.o arcos.o cartprint.o chainbuild...'
       write(iout,*)'GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES ...'
-      write(iout,*)'GAB: BIN = ../bin/unres_ifort_MPICH-restr-DFA_G...'
+      write(iout,*)'GAB: BIN = ../../../bin/unres/MD/unres_ifort_MP...'
       write(iout,*)'E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNR...'
       write(iout,*)'E0LL2Y: BIN = ../../../bin/unres/MD/unres_ifort...'
       write(iout,*)'++++ End of compile info ++++'

diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F
index 5707b4b..7cf67db 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -5989,7 +5989,8 @@ C AL 5/2/14 - Introduce list of restraints
          dij=dist(i,j)
          do k=1,constr_homology
            distance(k)=odl(k,ii)-dij
-           distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii)
+           distancek(k)=
+     &        0.5d0*waga_dist(iset)*distance(k)**2*sigma_odl(k,ii)
          enddo
          
          min_odl=minval(distancek)
@@ -6002,7 +6003,7 @@ C AL 5/2/14 - Introduce list of restraints
          odleg2=0.0d0
          do k=1,constr_homology
 c Nie wiem po co to liczycie jeszcze raz!
-c            odleg3=-waga_dist*((distance(i,j,k)**2)/ 
+c            odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ 
 c     &              (2*(sigma_odl(i,j,k))**2))
             godl(k)=dexp(-distancek(k)+min_odl)
             odleg2=odleg2+godl(k)
@@ -6023,8 +6024,8 @@ c Gradient
          sum_sgodl=0.0
          do k=1,constr_homology
 c            godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
-c     &           *waga_dist)+min_odl
-           sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
+c     &           *waga_dist(iset))+min_odl
+           sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist(iset)
            sum_sgodl=sum_sgodl+sgodl
 
 c            sgodl2=sgodl2+sgodl
@@ -6076,7 +6077,7 @@ c     &                                   -(6.28318-dih_diff(i,k))
 c          if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
 c     &                                   6.28318+dih_diff(i,k)
 
-          kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
+          kat3=-0.5d0*waga_angle(iset)*dih_diff(k)**2*sigma_dih(k,i)
           gdih(k)=dexp(kat3)
           kat2=kat2+gdih(k)
 c          write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
@@ -6099,7 +6100,7 @@ c ----------------------------------------------------------------------
         sum_gdih=kat2
         sum_sgdih=0.0
         do k=1,constr_homology
-          sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
+          sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle(iset)
           sum_sgdih=sum_sgdih+sgdih
         enddo
         grad_dih3=sum_sgdih/sum_gdih

diff --git a/source/unres/src_MD/geomout.F b/source/unres/src_MD/geomout.F
index a5c6f96..ac2c759 100644 (file)
--- a/source/unres/src_MD/geomout.F
+++ b/source/unres/src_MD/geomout.F
@@ -461,7 +461,7 @@ c-----------------------------------------------------------------
      &      (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair
      &             +21*nfrag_back
-        elseif(hremd.gt.0) then
+        elseif(hremd.gt.0.or.homol_nset.gt.1) then
            write(line2,'(i5)') iset
            format2="a005"
         else

diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F
index d21d3b9..18e2c2e 100644 (file)
--- a/source/unres/src_MD/readrtns.F
+++ b/source/unres/src_MD/readrtns.F
@@ -1166,6 +1166,8 @@ c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 
       if (constr_homology.gt.0) then
         call read_constr_homology
+      else
+        homol_nset=0
       endif
 
 
@@ -2412,6 +2414,7 @@ c-------------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
       include 'COMMON.MD'
+      include 'COMMON.CONTROL'
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath
       do i=1,2*nres
@@ -2420,7 +2423,7 @@ c-------------------------------------------------------------------------------
       do i=1,2*nres
          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
       enddo
-      if(usampl) then
+      if(usampl.or.homol_nset.gt.1) then
              read (irest2,*) iset
       endif
       close(irest2)
@@ -2637,12 +2640,26 @@ c-------------------------------------------------------------------------------
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
+      character*3200 controlcard1
       integer ki, i, j, k, l
       logical lprn /.true./
 
       call card_concat(controlcard)
-      call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
-      call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
+      call readi(controlcard,"HOMOL_NSET",homol_nset,1)       
+      if (homol_nset.gt.1)then
+         call card_concat(controlcard)
+         read(controlcard,*) (waga_dist(i),i=1,homol_nset) 
+         call card_concat(controlcard)
+         read(controlcard,*) (waga_angle(i),i=1,homol_nset) 
+         write(iout,*) "iset distance_weight angle_weight"
+         do i=1,homol_nset
+           write(iout,*) i,waga_dist(i),waga_angle(i)
+         enddo
+      else
+       iset=1
+       call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
+       call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
+      endif
 
       write (iout,*) "nnt",nnt," nct",nct
       call flush(iout)