From 5db6a4ce7482c069ab555ba9026247b8c3d69409 Mon Sep 17 00:00:00 2001 From: marsamos Date: Tue, 31 Aug 2010 03:22:54 +0000 Subject: [PATCH] path related routines in pwtools moved to NEB dir git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@7043 c92efa57-630b-4861-b058-cf58834340f0 --- NEB/Makefile | 7 +- NEB/make.depend | 7 + NEB/path_int.f90 | 359 ++++++++++++++++++++++++++++++++++++++++++++++ NEB/path_int.sh | 191 ++++++++++++++++++++++++ NEB/path_merge.sh | 132 +++++++++++++++++ 5 files changed, 695 insertions(+), 1 deletion(-) create mode 100644 NEB/path_int.f90 create mode 100755 NEB/path_int.sh create mode 100755 NEB/path_merge.sh diff --git a/NEB/Makefile b/NEB/Makefile index 93dd8dbed..b9e420b38 100644 --- a/NEB/Makefile +++ b/NEB/Makefile @@ -31,13 +31,18 @@ EEOBJS=../EE/libee.a TLDEPS=bindir mods libs liblapack libblas pw eelib -all : tldeps neb.x +all : tldeps neb.x path_int.x neb.x : $(NEBOBJS) libneb.a $(LIBOBJS) $(PWOBJS) $(QEMODS) $(LD) $(LDFLAGS) -o $@ \ $(NEBOBJS) libneb.a $(PWOBJS) $(EEOBJS) $(QEMODS) $(LIBOBJS) $(LIBS) - ( cd ../bin; ln -fs ../NEB/$@ . ) +path_int.x : path_int.o $(PWOBJS) $(QEMODS) $(LIBOBJS) + $(LD) $(LDFLAGS) -o $@ \ + path_int.o $(PWOBJS) $(QEMODS) $(LIBOBJS) $(LIBS) + - ( cd ../bin ; ln -fs ../NEB/$@ . ) + libneb.a : $(NEBLIBS) $(AR) $(ARFLAGS) $@ $? $(RANLIB) $@ diff --git a/NEB/make.depend b/NEB/make.depend index 51ebc6580..150373ed4 100644 --- a/NEB/make.depend +++ b/NEB/make.depend @@ -81,6 +81,13 @@ path_base.o : path_io_routines.o path_base.o : path_opt_routines.o path_base.o : path_reparametrisation.o path_base.o : path_variables.o +path_int.o : ../Modules/basic_algebra_routines.o +path_int.o : ../Modules/cell_base.o +path_int.o : ../Modules/constants.o +path_int.o : ../Modules/io_files.o +path_int.o : ../Modules/kind.o +path_int.o : ../Modules/splinelib.o +path_int.o : path_formats.o path_io_routines.o : ../Modules/basic_algebra_routines.o path_io_routines.o : ../Modules/cell_base.o path_io_routines.o : ../Modules/constants.o diff --git a/NEB/path_int.f90 b/NEB/path_int.f90 new file mode 100644 index 000000000..6db99604b --- /dev/null +++ b/NEB/path_int.f90 @@ -0,0 +1,359 @@ +! +! Copyright (C) 2004 PWSCF group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +MODULE int_global_variables + ! + USE kinds, ONLY : DP + ! + IMPLICIT NONE + ! + INTEGER :: N, dim + INTEGER :: old_num_of_images, new_num_of_images + INTEGER :: first_image, last_image + REAL (DP), ALLOCATABLE :: old_pos(:,:), new_pos(:,:) + REAL (DP), ALLOCATABLE :: old_PES_gradient(:,:), new_PES_gradient(:,:) + INTEGER, ALLOCATABLE :: fix_atom(:) + REAL (DP), ALLOCATABLE :: old_V(:), new_V(:) + REAL (DP), ALLOCATABLE :: d_R(:) + REAL (DP), ALLOCATABLE :: a(:), b(:), c(:), d(:), F(:) + REAL (DP), ALLOCATABLE :: old_mesh(:), new_mesh(:) + REAL (DP), ALLOCATABLE :: tangent(:) + CHARACTER(LEN=256) :: old_restart_file, new_restart_file + ! +END MODULE int_global_variables +! +! +PROGRAM images_interpolator + ! + USE kinds, ONLY : DP + USE io_files, ONLY : iunrestart + USE constants, ONLY : eps16 + USE path_formats + USE basic_algebra_routines, ONLY : norm + USE cell_base, ONLY : at, alat + USE int_global_variables + USE splinelib + ! + IMPLICIT NONE + ! + INTEGER :: i, j, ia + INTEGER :: istep, nstep, suspended_image + INTEGER :: ierr + REAL (DP) :: R, delta_R, x + LOGICAL :: no_interpolation + LOGICAL, EXTERNAL :: matches + CHARACTER (LEN=20) :: cell_parameters + CHARACTER (LEN=256) :: input_line + + ! + ! + ! ... the input file is read + ! + READ(*,*) N + READ(*,*) old_num_of_images + READ(*,*) new_num_of_images + READ(*,*) first_image + READ(*,*) last_image + READ(*,'(A)') old_restart_file + READ(*,'(A)') new_restart_file + READ(*,*) alat + ! + READ( UNIT = *, FMT = '(A)', IOSTAT = ierr ) cell_parameters + ! + IF ( ierr < 0 ) THEN + ! + WRITE(*,'(T2,"read_input: the card CELL_PARAMETERS not found")') + STOP + ! + ELSE IF ( ierr > 0 ) THEN + ! + WRITE(*,'(T2,"read_input: an error occured reading CELL_PARAMETERS")') + STOP + ! + END IF + ! + IF ( .NOT. matches( "CELL_PARAMETERS" , TRIM( cell_parameters ) ) ) THEN + ! + WRITE(*,'(T2,"read_input: ",A," is not a valid card")') cell_parameters + STOP + ! + END IF + ! + READ( UNIT = *, FMT = *, IOSTAT = ierr ) at + ! + IF ( ierr < 0 ) THEN + ! + WRITE(*,'(T2,"read_input: lattice vectors not found")') + STOP + ! + ELSE IF ( ierr > 0 ) THEN + ! + WRITE(*,'(T2,"read_input: an error occured reading laccice vectors")') + WRITE(*, FMT = lattice_vectors ) at + STOP + ! + END IF + ! + dim = 3 * N + ! + ALLOCATE( old_pos( dim, old_num_of_images ) ) + ALLOCATE( old_PES_gradient( dim, old_num_of_images ) ) + ALLOCATE( old_V( old_num_of_images ) ) + ALLOCATE( new_pos( dim, new_num_of_images ) ) + ALLOCATE( new_PES_gradient( dim, new_num_of_images ) ) + ALLOCATE( new_V( new_num_of_images ) ) + ALLOCATE( fix_atom( dim ) ) + ALLOCATE( d_R( dim ) ) + ALLOCATE( tangent( dim ) ) + ALLOCATE( a( old_num_of_images - 1 ) ) + ALLOCATE( b( old_num_of_images - 1 ) ) + ALLOCATE( c( old_num_of_images - 1 ) ) + ALLOCATE( d( old_num_of_images - 1 ) ) + ALLOCATE( F( old_num_of_images ) ) + ALLOCATE( old_mesh( old_num_of_images ) ) + ALLOCATE( new_mesh( new_num_of_images ) ) + ! + ! ... the old restart file is read + ! + OPEN( UNIT = iunrestart, FILE = old_restart_file, STATUS = "OLD", & + ACTION = "READ" ) + ! + no_interpolation = .FALSE. + ! + READ( UNIT = iunrestart, FMT = * ) ! RESTART INFORMATION + READ( UNIT = iunrestart, FMT = * ) istep + READ( UNIT = iunrestart, FMT = * ) nstep + READ( UNIT = iunrestart, FMT = * ) suspended_image + READ( UNIT = iunrestart, FMT = * ) ! conv_elec + ! + ! ... read either "ENERGIES, POSITIONS AND GRADIENTS" + ! + repeat_loop: DO + ! + READ( UNIT = iunrestart, FMT = '(256A)', IOSTAT = ierr ) input_line + ! + IF ( matches( input_line, & + "ENERGIES, POSITIONS AND GRADIENTS" ) ) EXIT repeat_loop + ! + IF ( ierr /= 0 ) THEN + ! + WRITE( *, '(/,5X,"an error occured reading",A)' ) old_restart_file + ! + STOP + ! + END IF + ! + END DO repeat_loop + ! + READ( UNIT = iunrestart, FMT = * ) ! Image: + READ( UNIT = iunrestart, FMT = * ) old_V(1) + ! + ia = 0 + ! + DO j = 1, dim, 3 + ! + ia = ia + 1 + ! + READ( UNIT = iunrestart, FMT = * ) & + old_pos(j,1), & + old_pos((j+1),1), & + old_pos((j+2),1), & + old_PES_gradient(j,1), & + old_PES_gradient((j+1),1), & + old_PES_gradient((j+2),1), & + fix_atom(j), & + fix_atom((j+1)), & + fix_atom((j+2)) + ! + old_PES_gradient(:,1) = old_PES_gradient(:,1) * fix_atom + ! + END DO + ! + DO i = 2, old_num_of_images + ! + READ( UNIT = iunrestart, FMT = * ) ! Image: + READ( UNIT = iunrestart, FMT = * ) old_V(i) + ! + DO j = 1, dim, 3 + ! + READ( UNIT = iunrestart, FMT = * ) & + old_pos(j,i), & + old_pos((j+1),i), & + old_pos((j+2),i), & + old_PES_gradient(j,i), & + old_PES_gradient((j+1),i), & + old_PES_gradient((j+2),i) + ! + END DO + ! + old_PES_gradient(:,i) = old_PES_gradient(:,i) * fix_atom + ! + END DO + ! + CLOSE( UNIT = iunrestart ) + ! + F = 0.D0 + ! + DO i = 2, ( old_num_of_images - 1 ) + ! + ! ... tangent to the path ( normalized ) + ! + tangent(:) = 0.5D0 * ( ( old_pos(:,i+1) - old_pos(:,i) ) / & + norm( old_pos(:,i+1) - old_pos(:,i) ) + & + ( old_pos(:,i) - old_pos(:,i-1) ) / & + norm( old_pos(:,i) - old_pos(:,i-1) ) ) + ! + tangent = tangent / norm( tangent ) + ! + F(i) = DOT_PRODUCT( - old_PES_gradient(:,i) , tangent ) + ! + END DO + ! + old_mesh(1) = 0.D0 + ! + DO i = 1, ( old_num_of_images - 1 ) + ! + d_R = old_pos(:,(i+1)) - old_pos(:,i) + ! + R = norm( d_R ) + ! + old_mesh(i+1) = old_mesh(i) + R + ! + a(i) = 2.D0 * ( old_V(i) - old_V(i+1) ) / R**(3) - & + ( F(i) + F(i+1) ) / R**(2) + ! + b(i) = 3.D0 * ( old_V(i+1) - old_V(i) ) / R**(2) + & + ( 2.D0 * F(i) + F(i+1) ) / R + ! + c(i) = - F(i) + ! + d(i) = old_V(i) + ! + END DO + ! + i = first_image + ! + delta_R = ( old_mesh(last_image) - & + old_mesh(first_image) ) / DBLE( new_num_of_images - 1 ) + ! + DO j = 0, ( new_num_of_images - 1 ) + ! + R = old_mesh( first_image ) + DBLE(j) * delta_R + ! + new_mesh(j+1) = R + ! + check_index: DO + ! + IF ( ( R > old_mesh(i+1) ) .AND. & + ( i < ( old_num_of_images - 1 ) ) ) THEN + ! + i = i + 1 + ! + ELSE + ! + EXIT check_index + ! + END IF + ! + END DO check_index + ! + x = R - old_mesh( i ) + ! + new_V(j+1) = a(i)*(x**3) + b(i)*(x**2) + c(i)*x + d(i) + ! + END DO + ! + new_mesh = new_mesh / old_mesh(old_num_of_images) + old_mesh = old_mesh / old_mesh(old_num_of_images) + ! + CALL dosplineint( old_mesh , old_pos , new_mesh , new_pos ) + ! + new_PES_gradient = 0.D0 + ! + new_PES_gradient(:,1) = old_PES_gradient(:,1) + ! + new_PES_gradient(:,new_num_of_images) = old_PES_gradient(:,old_num_of_images) + ! + ! ... the new restart file is written + ! + OPEN( UNIT = iunrestart, FILE = new_restart_file, STATUS = "UNKNOWN", & + ACTION = "WRITE" ) + ! + WRITE( UNIT = iunrestart, FMT = '("RESTART INFORMATION")' ) + ! + ! ... by default istep and nstep are set to zero + ! + WRITE( UNIT = iunrestart, FMT = '(I4)' ) 0 + WRITE( UNIT = iunrestart, FMT = '(I4)' ) 0 + WRITE( UNIT = iunrestart, FMT = '(I4)' ) 0 + WRITE( UNIT = iunrestart, FMT = '(A4)' ) 'F' + ! + WRITE( UNIT = iunrestart, FMT = '("ENERGIES, POSITIONS AND GRADIENTS")' ) + ! + DO i = 1, new_num_of_images + ! + WRITE( UNIT = iunrestart, FMT = '("Image: ",I4)' ) i + WRITE( UNIT = iunrestart, FMT = energy ) new_V(i) + ! + ia = 0 + ! + DO j = 1, dim, 3 + ! + ia = ia + 1 + ! + IF ( i == 1 ) THEN + ! + WRITE( UNIT = iunrestart, FMT = restart_first ) & + new_pos(j,i), & + new_pos((j+1),i), & + new_pos((j+2),i), & + new_PES_gradient(j,i), & + new_PES_gradient((j+1),i), & + new_PES_gradient((j+2),i), & + fix_atom(j), & + fix_atom((j+1)), & + fix_atom((j+2)) + ! + ELSE + ! + WRITE( UNIT = iunrestart, FMT = restart_others ) & + new_pos(j,i), & + new_pos((j+1),i), & + new_pos((j+2),i), & + new_PES_gradient(j,i), & + new_PES_gradient((j+1),i), & + new_PES_gradient((j+2),i) + ! + END IF + ! + END DO + ! + END DO + ! + WRITE( UNIT = iunrestart, FMT = '("END")' ) + ! + CLOSE( UNIT = iunrestart ) + ! + IF ( ALLOCATED( old_pos ) ) DEALLOCATE( old_pos ) + IF ( ALLOCATED( old_PES_gradient ) ) DEALLOCATE( old_PES_gradient ) + IF ( ALLOCATED( old_V ) ) DEALLOCATE( old_V ) + IF ( ALLOCATED( new_pos ) ) DEALLOCATE( new_pos ) + IF ( ALLOCATED( new_PES_gradient ) ) DEALLOCATE( new_PES_gradient ) + IF ( ALLOCATED( new_V ) ) DEALLOCATE( new_V ) + IF ( ALLOCATED( fix_atom ) ) DEALLOCATE( fix_atom ) + IF ( ALLOCATED( d_R ) ) DEALLOCATE( d_R ) + IF ( ALLOCATED( tangent ) ) DEALLOCATE( tangent ) + IF ( ALLOCATED( a ) ) DEALLOCATE( a ) + IF ( ALLOCATED( b ) ) DEALLOCATE( b ) + IF ( ALLOCATED( c ) ) DEALLOCATE( c ) + IF ( ALLOCATED( d ) ) DEALLOCATE( d ) + IF ( ALLOCATED( F ) ) DEALLOCATE( F ) + IF ( ALLOCATED( old_mesh ) ) DEALLOCATE( old_mesh ) + IF ( ALLOCATED( new_mesh ) ) DEALLOCATE( new_mesh ) + ! +END PROGRAM images_interpolator diff --git a/NEB/path_int.sh b/NEB/path_int.sh new file mode 100755 index 000000000..32a6abbed --- /dev/null +++ b/NEB/path_int.sh @@ -0,0 +1,191 @@ +#!/bin/bash --noprofile +# +################################################################################ +## Copyright (C) 2004 Carlo Sbraccia. ## +## This file is distributed under the terms ## +## of the GNU General Public License. ## +## See http://www.gnu.org/copyleft/gpl.txt . ## +## ## +## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ## +## EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF ## +## MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND ## +## NONINFRINGEMENT. IN NO EVENT SHALL CARLO SBRACCIA BE LIABLE FOR ANY ## +## CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, ## +## TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE ## +## SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ## +################################################################################ +# make sure there is no locale setting creating unneeded differences. +LC_ALL=C +export LC_ALL +# +################################################################################ +## Set these variables according to your needs ## +################################################################################ +# +# ... root directory of the PWscf package +# +ROOT_DIR="" +# +# ... old and new restart file name +# ... ( usually old_prefix.neb and new_prefix.neb ) +# +old_restart_file="" +new_restart_file="" +# +# ... number of images of the old and of the new path +# +old_num_of_images= +new_num_of_images= +# +# ... interpolation is performed between the first_image and the last_image +# ... of the old path +# +first_image= +last_image= +# +# ... the number of atoms +# +nat= +# +# ... lattice parameter ( celldm(1) in the pw input file ) +# +alat=1.D0 +# +# ... interpolation is possible only in the ibrav=0 case +# ... this card is the same used in the pw input file +# +cat > CELL_PARAMETERS << EOF +CELL_PARAMETERS + 0.00000 0.00000 0.00000 + 0.00000 0.00000 0.00000 + 0.00000 0.00000 0.00000 +EOF +# +# ... optional informations: +# ... needed to gereate visualization files ( xyz and axsf formats ) +# +# ... put a symbol for each atomic specie that compose the system +# +list_of_atoms="" +# +# ... number of atoms of each specie +# ... (the sum of all N[i] must be equal to nat) +# +N[1]= +# +################################################################################ +################################################################################ +###### DO NOT MODIFY THE SCRIPT UNDER THESE LINES ###### +################################################################################ +################################################################################ +# +GAWK=$( which gawk ) +# +################################################################################ +## lattice vectors for gawk scripts ## +################################################################################ +a11=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==2) {printf "%12.8f", $1} }' ) +a12=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==2) {printf "%12.8f", $2} }' ) +a13=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==2) {printf "%12.8f", $3} }' ) +a21=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==3) {printf "%12.8f", $1} }' ) +a22=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==3) {printf "%12.8f", $2} }' ) +a23=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==3) {printf "%12.8f", $3} }' ) +a31=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==4) {printf "%12.8f", $1} }' ) +a32=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==4) {printf "%12.8f", $2} }' ) +a33=$( cat CELL_PARAMETERS | $GAWK '{ if (NR==4) {printf "%12.8f", $3} }' ) +# +################################################################################ +## the input file for the iterpolator code is generated ## +################################################################################ +# +if [ ! -f ${old_restart_file} ]; then + echo "Error: file ${old_restart_file} not fount"; exit +fi +# +cat > input << EOF +${nat} +${old_num_of_images} +${new_num_of_images} +${first_image} +${last_image} +${old_restart_file} +${new_restart_file} +${alat} +EOF +# +cat CELL_PARAMETERS | $GAWK '{ if ( NR == 1 ) { print }; if ( NR > 1 ) \ +{ printf " %12.8f %12.8f %12.8f\n", $1, $2, $3} }' >> input +# +# +$ROOT_DIR/bin/path_int.x < input +# +if [[ "${list_of_atoms}" != "" ]]; then + # + ############################################################################## + ## dynamical generation of the "from_restart_to_axfs.gawk" script ## + ############################################################################## + file="from_restart_to_axsf.gawk" +cat > ${file} << EOF +BEGIN{ + a_0 = 0.529177 ; count = -10000; + printf " ANIMSTEPS %3i \n", ${new_num_of_images} ; + printf " CRYSTAL \n" ; + printf " PRIMVEC \n" ; + printf " %16.10f %16.10f %16.10f \n", ${a11}*a_0, ${a12}*a_0, ${a13}*a_0; + printf " %16.10f %16.10f %16.10f \n", ${a21}*a_0, ${a22}*a_0, ${a23}*a_0; + printf " %16.10f %16.10f %16.10f \n", ${a31}*a_0, ${a32}*a_0, ${a33}*a_0; +} +{ +if ( \$0 == "RESTART INFORMATIONS" ) { next; next; next } +if ( \$0 == "ENERGY, POSITIONS AND GRADIENTS" ) { next } +if ( \$0 == "VELOCITIES" ) { exit } +if ( \$1 == "Image:" ) { + count = -1 ; + printf " PRIMCOORD %3i \n", \$2 ; + printf "%4i 1 \n", ${nat} ; + } +else { + count++; +EOF + # + ref1=0 + ref2=0 + # + index=0 + # + for atom in ${list_of_atoms}; do + # + index=$(( ${index} + 1 )) + # + ref1=$(( ${ref2} + 1 )) + ref2=$(( ${ref2} + ${N[${index}]} )) + # + echo " if ( count >= ${ref1} && count <= ${ref2} ) { " >> ${file} + echo " printf \"${atom} \"; " >> ${file} + echo " printf \" %16.10f \", \$1 * a_0 ; " >> ${file} + echo " printf \" %16.10f \", \$2 * a_0 ; " >> ${file} + echo " printf \" %16.10f \", \$3 * a_0 ; " >> ${file} + echo " printf \" %16.10f \", \$4 / a_0 ; " >> ${file} + echo " printf \" %16.10f \", \$5 / a_0 ; " >> ${file} + echo " printf \" %16.10f \", \$6 / a_0 ; " >> ${file} + echo " printf \" \\n\"; " >> ${file} + echo " } " >> ${file} + # + done + # + echo " } " >> ${file} + echo "} " >> ${file} + # + ############################################################################## + # + $GAWK -f from_restart_to_axsf.gawk ${new_restart_file} > \ + ${new_restart_file}.axsf + # + rm -f from_restart_to_axsf.gawk + # +fi +# +echo "done" +# +rm -f input +rm -f CELL_PARAMETERS diff --git a/NEB/path_merge.sh b/NEB/path_merge.sh new file mode 100755 index 000000000..8730832b1 --- /dev/null +++ b/NEB/path_merge.sh @@ -0,0 +1,132 @@ +#!/bin/bash +################################################################################ +## Copyright (C) 2007 Guido Fratesi. ## +## This file is distributed under the terms ## +## of the GNU General Public License. ## +## See http://www.gnu.org/copyleft/gpl.txt . ## +## ## +## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ## +## EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF ## +## MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND ## +## NONINFRINGEMENT. IN NO EVENT SHALL GUIDO FRATESI BE LIABLE FOR ANY ## +## CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, ## +## TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE ## +## SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ## +################################################################################ +# # +# Suitable for PWSCF v. 4.0, might require revisions for future releases. # +# Execute with no arguments for help. # +# # +################################################################################ + +function showhelp () { +cat < /dev/stderr + + USAGE: path_merge.bash file0.path n0 m0 file1.path n1 m1 + + merges the two path files file0.path and file1.path by taking the + images n0 to m0 from file0.path and images n1 to m1 from file file1.path. + Merged file is written to standard output. + +EOF + +if [ "$err" ]; then cat << EOF > /dev/stderr + +* ERROR: $err + +EOF +fi +exit +} + +awk=`which awk` + +# +# ... Parse and check input parameters +# +inpf[0]=$1 # first path file (file0.path) +imlw[0]=$2 # from this image index (n0) +imup[0]=$3 # to this image index (m0) +inpf[1]=$4 # second path file (file1.path) +imlw[1]=$5 # from this image index (n1) +imup[1]=$6 # to this image index (m1) +# +if (($#==0)); then err=""; showhelp; fi +if (($#<6)); then err="Too few arguments!"; showhelp; fi +if (($#>6)); then err="Too many arguments!"; showhelp; fi +for i in 0 1; do + if [ ! -f ${inpf[$i]} ]; then + err="File n.$i \"${inpf[$i]}\" not found!" + showhelp + fi + # + # ... Extract number of atoms + # + nat[$i]=`$awk ' + ($1=="Image:")&&($2==1) {n=NR}; + ($1=="Image:")&&($2==2) {printf NR-n-2; exit}; + ' ${inpf[$i]}` + # + # ... Read number of images + # + nim[$i]=`$awk '/NUMBER OF IMAGES/ {getline; printf $1; exit}' ${inpf[$i]}` + # + if ((${imlw[$i]}<1)) || ((${imlw[$i]}>${nim[$i]})) || + ((${imup[$i]}<1)) || ((${imup[$i]}>${nim[$i]})) || + ((${imup[$i]}<${imlw[$i]})) ; then + err="Check the images requested from file n.$((i+1))" + showhelp + fi + # +done +# +if ((${nat[0]}!=${nat[1]})); then + err="The number of atoms in the two path files is not the same!" + showhelp +fi + + +# +# ... Write the header +# +cat <=nfix+1)&&(NR<=nfix+nat) {fix[NR-nfix]=sprintf("%3i%3i%3i",$7,$8,$9)}; + # + # ... Write images in the range requested + # + ($1=="Image:")&&($2>=imlw)&&($2<=imup) { n=NR+1; + printf("Image:%5i\n",iim); } + (NR==n) + (NR>=n+1)&&(NR<=n+nat) { for(i=1;i<=6;i++) printf("%20.12f",$i); + if(iim==1) printf("%9s",fix[NR-n]); + printf("\n"); }; + (NR==n+nat) { iim++; }; + ' ${inpf[$i]} + # + iim=$(($iim+${imup[$i]}-${imlw[$i]}+1)) + # +done +