program make_images
	USE IMAGE_MODULE 
        implicit none
!----------------------------------------------------------------------------------------------------------------------------------
!	The purpose of this program is to read  a starting and
!	final geometry, 
!	calculate a position vector that is the difference
!	between them, and use it to interpolate between the
!	starting and final geometries to generate a series
! 	of images suitable for use with the vasp nudged elastic
!	band program.
!
!	convention used for numbering the images:
!	00  		starting coordinates
!	01  		first image
!	02  		second image
!             ...
!	num_images 	last image
!	num_images + 1	product coordinates
!
!
!----------------------------------------------------------------------------------------------------------------------------------
!	Declarations:
!----------------------------------------------------------------------------------------------------------------------------------
	integer :: num_images
	integer :: counter  
	type (model) :: mod 
	type (image) :: reactant, product 
	type (vector)  :: vec
	character (len=20) :: r_file, p_file, this_file
	print *, "assumes non-neg scaling factor" 
	print *, "assumes positions given in direct lattice"
	print *, "The number of images must be 2^n"
	print *, "How many images do you want to make?"
	read(unit=*,fmt=*) num_images
	print *, "How many species are present in the model?"
	read(unit=*,fmt=*) mod%num_spec
	print *, "What is the total number of ions present?"
	read(unit=*,fmt=*) mod%t_num_ions
	print *, "what is the name of the file for the reactant?"
	read(unit=*,fmt=*) r_file 
        print *, "what is the name of the file for the product?"
	read(unit=*,fmt=*) p_file 
        call read_POSCAR(mod, reactant, r_file) 
       	call read_POSCAR(mod, product, p_file)
        call make_difference_vector(mod, reactant, product, vec)
        call divide_difference_vector(mod, vec, num_images)
	counter = 0
	do while (counter <= (num_images+1))
	call write_POSCAR(counter, mod, reactant, vec)
	  counter = counter + 1
	end do
	end program make_images
!----------------------------------------------------------------------------------------------------------------------------------
	subroutine read_POSCAR(mod, this_im, this_file)
	USE IMAGE_MODULE
	implicit none
	type (model) :: mod 
        type (image) :: this_im
	integer :: a,b
	character (len=20) :: this_file
	print *, "********************************************"
	print *, "inside read_POSCAR"
	open(unit=15,file=this_file,action="read", status="old")
	read(unit=15,fmt=*) this_im%name 
	write(unit=*,fmt=*)"name is ", this_im%name
	read(unit=15,fmt=*) this_im%scale 
	write(unit=*,fmt=*)"scale is ", this_im%scale
	if (this_im%scale.le.0) then
	write(unit=*,fmt=*)"scaling factor is zero or negative"
	end if
	read(unit=15,fmt=*)((this_im%lattvec(a,b),b=1,3),a=1,3)
	do a=1,3
	write(unit=*,fmt=3)"latt_vec",(this_im%lattvec(a,b),b=1,3)
	end do
 3	format(A10,2x,3(1x,f18.15))
 4	format(A10,2x,3(1x,f18.15),A10,2x,3(1x,A2))
	read(unit=15,fmt=*)(this_im%num_atoms_of_species(a),a=1,mod%num_spec)
	write(unit=*,fmt=5)(this_im%num_atoms_of_species(a),a,a=1,mod%num_spec)
 5	format(1x, I2,' atoms of species',I2)
	read(unit=15,fmt=*) this_im%seldym
	read(unit=15,fmt=*) this_im%direct
	write(unit=*,fmt=*)'D if coords are direct: ', this_im%direct
	do a=1,mod%t_num_ions
	read(unit=15,fmt=*)(this_im%coords(a,b),b=1,3),(this_im%dyn(a,b),b=1,3)
	end do
	do a=1,mod%t_num_ions
	write(unit=*,fmt=4)"coord",(this_im%coords(a,b),b=1,3),"dynamics",(this_im%dyn(a,b),b=1,3)
	end do
	close(unit=15)
	return
	end 	subroutine read_POSCAR
!----------------------------------------------------------------------------------------------------------------------------------
	subroutine make_difference_vector(mod, reactant, product, vec)
	USE IMAGE_MODULE
	implicit none
	type (model) :: mod 
	type (image) :: reactant, product
	type (vector) :: vec
	integer :: a,b
	print *, "********************************************"
        print *, "inside make_difference_vector" 
	vec%direction = product%coords - reactant%coords
	do a=1,mod%t_num_ions
	write(unit=*,fmt=3)"direction",(vec%direction(a,b),b=1,3)
  3	format(A10,2x,3(1x,f18.16))
	end do
	return
	end 	subroutine make_difference_vector
!----------------------------------------------------------------------------------------------------------------------------------
	subroutine divide_difference_vector(mod, vec, num_images) 
	USE IMAGE_MODULE
	implicit none
	type (model) :: mod 
	type (vector) :: vec
	integer :: num_images
        print *, "********************************************"
        print *, "inside divide_difference_vector"
	vec%fraction = vec%direction /(num_images + 1)
	return
	end 	subroutine divide_difference_vector
!----------------------------------------------------------------------------------------------------------------------------------
	subroutine write_POSCAR(counter, mod, reactant, vec)
	USE IMAGE_MODULE
	implicit none
	type (model) :: mod 
        type (image) :: reactant, this_im 
        type (vector) :: vec
	integer, intent(in) :: counter 
	integer  :: a,b,num 
        print *, "********************************************"
        print *, "inside write_POSCAR"
	print *, counter 
	num = counter+10
	this_im%name                 = reactant%name
	this_im%direct               = reactant%direct
	this_im%seldym	             = reactant%seldym	
	this_im%lattvec              = reactant%lattvec
	this_im%scale                = reactant%scale
	this_im%num_atoms_of_species = reactant%num_atoms_of_species
	this_im%sum                  = reactant%sum
	this_im%dyn		     = reactant%dyn	
	this_im%coords = reactant%coords + (counter*vec%fraction)
        do a=1,mod%t_num_ions
	write(unit=*,fmt=4)"coord",(this_im%coords(a,b),b=1,3),"dynamics",(this_im%dyn(a,b),b=1,3)
 4      format(A10,2x,3(1x,f18.16),A10,2x,3(1x,A2))
        end do
	open(unit=num, action="write", status="unknown")
	write(unit=num, fmt=7) this_im%name, counter
 7 	format(A10,2x,'image number',I2)
	write(unit=num, fmt=3) this_im%scale
 3	format(f18.15)
	write(unit=*, fmt=*) this_im%scale
	write(unit=num, fmt=5) ((this_im%lattvec(a,b),b=1,3),a=1,3)
        write(unit=*, fmt=5) ((this_im%lattvec(a,b),b=1,3),a=1,3)
        write(unit=*, fmt=*) ((this_im%lattvec(a,b),b=1,3),a=1,3)
 5      format(3(1x,f18.15))
	write(unit=num, fmt=8) (this_im%num_atoms_of_species(a),a=1,mod%num_spec) 
	write(unit=num, fmt=6) this_im%seldym
	write(unit=num, fmt=6) this_im%direct
 6      format(A6)
	do a=1,mod%t_num_ions
	write(unit=num, fmt=9)(this_im%coords(a,b),b=1,3),(this_im%dyn(a,b),b=1,3)
	end do 
 8 	format(20(1x,I2))
 9      format(3(1x,f18.16),3(1x,A2))

	close(unit=20)
	return
	end 	subroutine write_POSCAR 
!----------------------------------------------------------------------------------------------------------------------------------
        MODULE IMAGE_MODULE

	type model
	integer :: num_spec
	integer :: t_num_ions 
        end type model

        type image
	character (len=20) ::  name
        character (len=6)  ::  direct, seldym
        real  :: scale
        integer :: a,b
        real,dimension(1:3,1:3) :: lattvec
        integer,dimension(1:2)  :: num_atoms_of_species
 	real,dimension(1:48,1:3)     :: sum
        real,dimension(1:48,1:3)     :: coords
	real,dimension(1:48,1:3)     :: coords_vec
	real,dimension(1:48,1:3)     :: coords_scale 
	real,dimension(1:3)         :: sum_pos
	character(len=2),dimension(1:48,1:3)           :: dyn
        end type image

	type vector 
	real,dimension(1:48,1:3)     :: direction
	real,dimension(1:48,1:3)     :: fraction
	end type vector 
        end MODULE IMAGE_MODULE