From 624886bbe93f5929cf6271438dac33374531eb27 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Wed, 25 Sep 2019 20:11:10 -0400 Subject: [PATCH 1/8] Initial commit with working atom and lattice type parsing --- README.md | 96 ++++++++++++++++++++++++++++++++++- example.in | 7 +++ src/Makefile | 18 +++++++ src/elements.f90 | 32 ++++++++++++ src/lattice.f90 | 95 ++++++++++++++++++++++++++++++++++ src/main.f90 | 52 +++++++++++++++++++ src/precision_comm_module.f90 | 13 +++++ src/region.f90 | 14 +++++ src/subroutines.f90 | 24 +++++++++ 9 files changed, 350 insertions(+), 1 deletion(-) create mode 100644 example.in create mode 100644 src/Makefile create mode 100644 src/elements.f90 create mode 100644 src/lattice.f90 create mode 100644 src/main.f90 create mode 100644 src/precision_comm_module.f90 create mode 100644 src/region.f90 create mode 100644 src/subroutines.f90 diff --git a/README.md b/README.md index 9cc824a..4b5280c 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,96 @@ # CAC_Model_Builder -This is a tool for building models in CAC +This is a tool for building models in CAC. Commands and usage options are below. + + + +## Flow of commands + +```flow +op1=>operation: Define atom types and lattices to be used +op2=>operation: Define regions and build +op3=>operation: Define modifiers +op4=>operation: Output data files +op1->op2->op3->op4 +``` + + + +## Command syntax + +### Atom types command + +``` +atom_types num_atoms {name mass} +``` + +The parameters for the atoms command are: + +`num_atoms` - number of atom types defined for this model building session + +`{}` - indicate that the contents must be repeated `num_atoms` times. + +`name` - Elemental name of atom + +`mass` - mass of the atom + +This command should only be called once, defining all atoms in one go. The atom types will then be defined in numeric order with the first atom defined being type one and the last one being type `num_atoms`. + +### Lattice command + +``` +lattice id lattice_type lattice_parameter [type atom_type] [basis num_basis_atoms {type posx posy posz}] +``` + +The parameters for the lattice command are: + +`id` - User defined id for this lattice type + +`lattice_type` - One of predefined lattice types which specifies the element type used. Current accepted options are: `FCC` + +`type` - Optional keyword which defines the atom type used for the lattice. This is used in place of basis if atoms are at lattice positions in these elements. + +`atom_type` - The atom type which corresponds to the atoms at the lattice positions of the current element + +`basis` - Optional keyword which is used in order to define the basis atoms instead of using the default definition. If basis is not included the following commands also are not included. + +`num_basis_atoms` are the number of basis atoms in this element. + +`{}` - indicate that the contents must be repeated `num_basis_atoms` times. + +`type` - the atom type of the atom. + +`posx posy posz` - The position of the basis atom relative to the lattice point at (0,0,0) + +**Either type or basis keywords must be present in the lattice command, both cannot be used.** + +## Region Command + +``` +region id lattice_id element_size units lenx leny lenz [zigzag] [origin x y z] [cat region_id dim [nomatch]] [orient [hkl] [hkl] [hkl]] +``` + +`id` - User defined id for this region + +`lattice_id` - The lattice type for this region + +`element_size` - The element size used for this region defined as the number of atoms per element edge. An element size of 2 means that this region is at full atomistic resolution. + +`units` - Either `lattice` or `box` which adjusts how the length values are calculated. Units `lattice` means the region will consist of `len` number of elements for every dim. Units `box` are defined in angstroms. + +`lenx leny lenz` - The lengths of the box in each dimension in the user defined units + +`zigzag` - Optional keyword which specifies if regions built with elements should have filled in boundaries (using atoms). If zigzag isn't present then the regions are built with filled in boundaries by default + +`origin x y z` - Optional keyword which specifies the origin of the current region in angstroms. The region boundaries are then (x, x+lenx), (y, y+leny), (z,z+lenz). + +`cat region_id outregionid dim [nomatch] ` - Optional keyword which stacks the current region on the face of another region defined by `dim`. `region_id` is the id of the region which is already build. `outregionid` is the user defined id of the combined stacked region which can be used with further merge commands. Default behavior is to expand the smallest region to match the larger one, using the optional keyword `nomatch` preserves the original regions and does not attempt to match the boundaries. + +`orient [hkl] [hkl] [hkl]` simply orients the unit cell of this region. This defaults to [100] [010] [001] + +### Write command + +``` +write file_name +``` + +Self explanatory. \ No newline at end of file diff --git a/example.in b/example.in new file mode 100644 index 0000000..f9cf796 --- /dev/null +++ b/example.in @@ -0,0 +1,7 @@ +#This is an example input script for the CAC model builder + +atom_types 1 Cu 63.546 +lattice 1 fcc 3.615 type 1 + +#region 1 1 2 lattice 20 20 20 +#write atoms.xyz \ No newline at end of file diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..9acc4e5 --- /dev/null +++ b/src/Makefile @@ -0,0 +1,18 @@ +CC=ifort +FC=ifort +FFLAGS=-c -mcmodel=large -debug -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone + +OBJECTS= main.o elements.o lattice.o subroutines.o precision_comm_module.o + +.SUFFIXES: +.SUFFIXES: .c .f .f90 .F90 .o + +builder: $(OBJECTS) + $(FC) $(OBJECTS) -o $@ + +.f90.o: + $(FC) $(FFLAGS) $< + +main.o lattice.o elements.o region.o subroutines.o : precision_comm_module.o +lattice.o : subroutines.o +main.o : elements.o lattice.o region.o \ No newline at end of file diff --git a/src/elements.f90 b/src/elements.f90 new file mode 100644 index 0000000..1fbacf6 --- /dev/null +++ b/src/elements.f90 @@ -0,0 +1,32 @@ +module elements + + use precision_comm_module + + implicit none + + !This is the data structure which is used to represent the CAC elements + type element + integer :: tag = 0 !Element tag (used to keep track of id's + integer :: type = 0 !Lattice type of the element + integer :: size = 0 !Element size + !Nodal position array below only works for wedge or fcc elements + real(kind=wp) :: r_node(3,8) + end type + + !Finite element array + type(element), allocatable :: element_array(:) + + integer :: ele_num + + !Data structure used to represent atoms + type atom + integer :: tag = 0 + integer :: type = 0 + real(kind =wp) :: r + end type + + type(atom), allocatable :: atoms(:) + + integer :: atom_num + +end module elements \ No newline at end of file diff --git a/src/lattice.f90 b/src/lattice.f90 new file mode 100644 index 0000000..4783a2b --- /dev/null +++ b/src/lattice.f90 @@ -0,0 +1,95 @@ +module lattice + + use precision_comm_module + use subroutines + + implicit none + + integer :: atom_types + + !Atom type variables + character(len=2), dimension(10) :: atom_names + real(kind=wp), dimension(10) :: atom_masses + + !Lattice_type variables + integer :: lat_num + character(len=10), dimension(10) :: lattice_id, lattice_type + real(kind=wp), dimension(10) :: lapa + integer(kind=wp), dimension(10) :: basis_atom_num + integer(kind=wp), dimension(10,10) :: basis_type + real(kind=wp), dimension(3,10,10) ::basis_pos + + !Unit Cell variables + real(kind = wp) :: fcc_cell(3,8), fcc_mat(3,3) + + public + contains + + subroutine lattice_init + + !Initialize needed variables + lat_num=0 + basis_atom_num(:) = 0 + + !Initialize finite element cells to be used + + !First initialize the primitive fcc cell + fcc_cell = reshape((/ 0.0_wp, 0.0_wp, 0.0_wp, & + 0.5_wp, 0.5_wp, 0.0_wp, & + 0.5_wp, 1.0_wp, 0.5_wp, & + 0.0_wp, 0.5_wp, 0.5_wp, & + 0.5_wp, 0.0_wp, 0.5_wp, & + 1.0_wp, 0.5_wp, 0.5_wp, & + 1.0_wp, 1.0_wp, 1.0_wp, & + 0.5_wp, 0.5_wp, 1.0_wp /), & + shape(fcc_cell)) + + fcc_mat = reshape((/ 0.5_wp, 0.5_wp, 0.0_wp, & + 0.5_wp, 0.5_wp, 0.5_wp, & + 0.5_wp, 0.0_wp, 0.5_wp /), & + shape(fcc_mat)) + end subroutine lattice_init + !This subroutine defines the atom type arrays + subroutine atom_type_parse(line) + + character(len=100), intent(in) :: line + character(len=100) :: errorloc + integer :: ia, error + character(len=20) :: label + + read(line, *, iostat=error) label, atom_types, (atom_names(ia), atom_masses(ia), ia=1, atom_types) + errorloc="lattice:22" + call read_error_check(error,errorloc) + + end subroutine atom_type_parse + + !This subroutine defines the lattice types and the unit cells for the lattice types + subroutine lattice_parse(line) + character(len=100), intent(in) :: line + integer :: ia, error + character(len=20) :: label, kw + character(len=100) :: errorloc + + lat_num = lat_num + 1 + read(line, *, iostat=error) label, lattice_id(lat_num), lattice_type(lat_num), lapa(lat_num), kw + errorloc="lattice:77" + call read_error_check(error, errorloc) + + select case(kw) + case("type") + read(line(scan(line, "type"):), *, iostat=error) label, basis_type(1,1) + errorloc="lattice:56" + call read_error_check(error,errorloc) + case("basis") + read(line(scan(line, "basis"):), *, iostat=error) label, basis_atom_num(lat_num), (basis_type(ia, lat_num) ,& + basis_pos(1:3,ia,lat_num), ia = 1, basis_atom_num(lat_num)) + errorloc="lattice:59" + call read_error_check(error,errorloc) + case default + print *, "Keyword ", kw, " is not accepted in the lattice command" + stop "Exit with error" + end select + + + end subroutine lattice_parse +end module lattice \ No newline at end of file diff --git a/src/main.f90 b/src/main.f90 new file mode 100644 index 0000000..323f996 --- /dev/null +++ b/src/main.f90 @@ -0,0 +1,52 @@ +program main + + use precision_comm_module + use elements + use lattice + use region + + integer :: iosline, iospara + logical :: flags(4) + character(len=100) :: line, command, errorloc + + iosline = 0 + iospara = 0 + flags(:) = .false. + + call lattice_init + + !Main command loop + do while (iosline == 0) + + read(*, '(a)', iostat=iosline) line + errorloc="read_input:line" + call read_error_check(iosline, errorloc) + + !Check for comment character (#) + if ((scan(line, '#')/= 1).and.(line/='')) then + read(line, *, iostat = iospara) command + errorloc="read_input:command" + call read_error_check(iosline, errorloc) + + select case(command) + case('atom_types') + call atom_type_parse(line) + flags(1) = .true. + case('lattice') + if(flags(1).eqv..false.) then + print *, "Please define atom types before defining lattice types" + stop 3 + end if + call lattice_parse(line) + flags(2) =.true. + ! case('region') + ! call build_region(line) + ! case('write') + ! call write_parse(line) + case default + print *, "The command ", trim(command), " is not currently accepted",& + " please check input script and try again." + end select + end if + end do +end program main \ No newline at end of file diff --git a/src/precision_comm_module.f90 b/src/precision_comm_module.f90 new file mode 100644 index 0000000..7532233 --- /dev/null +++ b/src/precision_comm_module.f90 @@ -0,0 +1,13 @@ + module precision_comm_module + + implicit none + + integer, parameter :: & + dp = selected_real_kind(15, 307), & ! double real + qp = selected_real_kind(33, 4931), & ! quadrupole real + wp = dp + + integer, save :: & + mpi_wp + + end module precision_comm_module diff --git a/src/region.f90 b/src/region.f90 new file mode 100644 index 0000000..f1a9d33 --- /dev/null +++ b/src/region.f90 @@ -0,0 +1,14 @@ +module region + use precision_comm_module + + implicit none + + public + contains + + subroutine build_region(line) + + character(len=100), intent(in) :: line + + end subroutine build_region +end module region \ No newline at end of file diff --git a/src/subroutines.f90 b/src/subroutines.f90 new file mode 100644 index 0000000..71aa475 --- /dev/null +++ b/src/subroutines.f90 @@ -0,0 +1,24 @@ +module subroutines + + use precision_comm_module + + implicit none + + public + contains + + + !This subroutine is just used to break the code and exit on an error + subroutine read_error_check(para, loc) + + integer, intent(in) :: para + character(len=100), intent(in) :: loc + + if (para > 0) then + print *, "Read error in ", trim(loc), " because of ", para + stop "Exit with error" + end if + end subroutine + + +end module subroutines \ No newline at end of file From 8217e8b51c61f5fa8ae7493a3e017cf3cd53887f Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 25 Nov 2019 18:19:25 -0500 Subject: [PATCH 2/8] Initial commit with file writing and create mode working for atoms --- README.md | 102 ++- example.in | 7 - src/Makefile | 32 +- src/atoms.f90 | 1167 +++++++++++++++++++++++++++++++++ src/call_mode.f90 | 24 + src/elements.f90 | 197 +++++- src/functions.f90 | 235 +++++++ src/io.f90 | 177 +++++ src/lattice.f90 | 95 --- src/main.f90 | 71 +- src/mode_create.f90 | 432 ++++++++++++ src/parameters.f90 | 8 + src/precision_comm_module.f90 | 13 - src/region.f90 | 14 - src/subroutines.f90 | 132 +++- 15 files changed, 2450 insertions(+), 256 deletions(-) delete mode 100644 example.in create mode 100644 src/atoms.f90 create mode 100644 src/call_mode.f90 create mode 100644 src/functions.f90 create mode 100644 src/io.f90 delete mode 100644 src/lattice.f90 create mode 100644 src/mode_create.f90 create mode 100644 src/parameters.f90 delete mode 100644 src/precision_comm_module.f90 delete mode 100644 src/region.f90 diff --git a/README.md b/README.md index 4b5280c..9fa35ee 100644 --- a/README.md +++ b/README.md @@ -1,96 +1,94 @@ # CAC_Model_Builder -This is a tool for building models in CAC. Commands and usage options are below. +This is a tool for building models in CAC. Commands and usage options are below. This code is intended to follow the atomsk code fairly closely. +## Modes +The modes follow similarly to the modes you find when using atomsk. The modes will be listed below alongside their syntax and other usage instructions. As a note, if a mode is being used then it has to come first. -## Flow of commands +### Mode Create -```flow -op1=>operation: Define atom types and lattices to be used -op2=>operation: Define regions and build -op3=>operation: Define modifiers -op4=>operation: Output data files -op1->op2->op3->op4 +``` +cacmb --create name element_type lattice_parameter esize ``` +Mode create has the following parameters: +`name` - User defined name that either defines the atom type or the lattice type if using the basis option -## Command syntax +`element_type` - Specifies which element type to use, this dictates the crystal being build. Current acceptable options for element_type are: -### Atom types command +* FCC - Uses the Rhombohedral primitive fcc unit cell as the finite element. -``` -atom_types num_atoms {name mass} -``` - -The parameters for the atoms command are: - -`num_atoms` - number of atom types defined for this model building session +`lattice_parameter` - The lattice parameter for the crystal structure. -`{}` - indicate that the contents must be repeated `num_atoms` times. +`esize` - Number of atoms per edge of the finite element. A value of 2 signifies full atomistic resolution and is the lowest number acceptable. -`name` - Elemental name of atom - -`mass` - mass of the atom - -This command should only be called once, defining all atoms in one go. The atom types will then be defined in numeric order with the first atom defined being type one and the last one being type `num_atoms`. - -### Lattice command +**Example** ``` -lattice id lattice_type lattice_parameter [type atom_type] [basis num_basis_atoms {type posx posy posz}] +cacmb --create Cu fcc 3.615 11 ``` -The parameters for the lattice command are: - -`id` - User defined id for this lattice type +Creates a copper element with a lattice parameter of 3.615 with 11 atoms per side -`lattice_type` - One of predefined lattice types which specifies the element type used. Current accepted options are: `FCC` +#### Optional keywords -`type` - Optional keyword which defines the atom type used for the lattice. This is used in place of basis if atoms are at lattice positions in these elements. +**Orient** -`atom_type` - The atom type which corresponds to the atoms at the lattice positions of the current element +``` orient [hkl] [hkl] [hkl] +orient [hkl] [hkl] [hkl] +``` -`basis` - Optional keyword which is used in order to define the basis atoms instead of using the default definition. If basis is not included the following commands also are not included. +Default orientation is `[100] [010] [001]`. If this keyword is present then the user must provide the orientation matrix in form `[hkl] [hkl] [hkl]`. -`num_basis_atoms` are the number of basis atoms in this element. +*Example:* `orient [-112] [110] [-11-1]` -`{}` - indicate that the contents must be repeated `num_basis_atoms` times. +**Basis** -`type` - the atom type of the atom. +``` +basis num atom_name x y z +``` -`posx posy posz` - The position of the basis atom relative to the lattice point at (0,0,0) +Default basis has `atom_name = name` with position (0,0,0). If used then the `atom_name x y z` must be include `num` times. -**Either type or basis keywords must be present in the lattice command, both cannot be used.** +*Example:* `basis 2 Mg 0 0 0 Mg 0.5 0.288675 0.81647` -## Region Command +**Duplicate** ``` -region id lattice_id element_size units lenx leny lenz [zigzag] [origin x y z] [cat region_id dim [nomatch]] [orient [hkl] [hkl] [hkl]] +duplicate numx numy numz ``` -`id` - User defined id for this region +Default duplicate is `1 1 1`. This is used to replicate the element along each dimensions. This cannot be used if the keyword dimensions is included. By default jagged edges along boundaries are filled if duplicate is greater than `1 1 1`. -`lattice_id` - The lattice type for this region +*Example:* `duplicate 10 10 10` -`element_size` - The element size used for this region defined as the number of atoms per element edge. An element size of 2 means that this region is at full atomistic resolution. +**Dimensions** + +``` +dimensions dimx dimy dimz +``` -`units` - Either `lattice` or `box` which adjusts how the length values are calculated. Units `lattice` means the region will consist of `len` number of elements for every dim. Units `box` are defined in angstroms. +There is no default dimensions as duplicate is the default option. This command assigns a box with user-assigned dimensions and fills it with the desired element. By default atoms fill in the jagged edges at the boundaries if the dimensions command is included. -`lenx leny lenz` - The lengths of the box in each dimension in the user defined units +Example: `dimensions 100 100 100` -`zigzag` - Optional keyword which specifies if regions built with elements should have filled in boundaries (using atoms). If zigzag isn't present then the regions are built with filled in boundaries by default +**ZigZag** -`origin x y z` - Optional keyword which specifies the origin of the current region in angstroms. The region boundaries are then (x, x+lenx), (y, y+leny), (z,z+lenz). +``` +zigzag boolx booly boolz +``` -`cat region_id outregionid dim [nomatch] ` - Optional keyword which stacks the current region on the face of another region defined by `dim`. `region_id` is the id of the region which is already build. `outregionid` is the user defined id of the combined stacked region which can be used with further merge commands. Default behavior is to expand the smallest region to match the larger one, using the optional keyword `nomatch` preserves the original regions and does not attempt to match the boundaries. +Default zigzag is `f f f`. This command specifies whether a boundary should be left jagged (i.e. in essence not filled in). If `boolx` is `t` than the x dimension is left jagged and if it is `f` then the x dimension is filled. -`orient [hkl] [hkl] [hkl]` simply orients the unit cell of this region. This defaults to [100] [010] [001] +*Example:* `zigzag t f t` gives a box with jagged edges in the x and z and filled edges in the y. -### Write command +**Origin** ``` -write file_name +origin x y z ``` -Self explanatory. \ No newline at end of file +Default origin is `0 0 0`. This command just sets the origin for where the simulation cell starts building. + +*Example:* `origin 10 0 1` \ No newline at end of file diff --git a/example.in b/example.in deleted file mode 100644 index f9cf796..0000000 --- a/example.in +++ /dev/null @@ -1,7 +0,0 @@ -#This is an example input script for the CAC model builder - -atom_types 1 Cu 63.546 -lattice 1 fcc 3.615 type 1 - -#region 1 1 2 lattice 20 20 20 -#write atoms.xyz \ No newline at end of file diff --git a/src/Makefile b/src/Makefile index 9acc4e5..c83df06 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,18 +1,34 @@ -CC=ifort FC=ifort -FFLAGS=-c -mcmodel=large -debug -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone - -OBJECTS= main.o elements.o lattice.o subroutines.o precision_comm_module.o +#FFLAGS=-c -mcmodel=large -debug -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone +FFLAGS=-c -mcmodel=large -Ofast +MODES=mode_create.o +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o build_subroutines.o $(MODES) .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o -builder: $(OBJECTS) +cacmb: $(OBJECTS) $(FC) $(OBJECTS) -o $@ .f90.o: $(FC) $(FFLAGS) $< -main.o lattice.o elements.o region.o subroutines.o : precision_comm_module.o -lattice.o : subroutines.o -main.o : elements.o lattice.o region.o \ No newline at end of file +.PHONY: clean +clean: + $(RM) cacmb *.o + +testfuncs: testfuncs.o functions.o subroutines.o + $(FC) testfuncs.o functions.o build_subroutines.o subroutines.o elements.o -o $@ + +.PHONY: cleantest +cleantest: + $(RM) testfuncs testfuncs.o + +$(OBJECTS) : parameters.o +atoms.o subroutines.o testfuncs.o : functions.o +main.o io.o build_subroutines.o: elements.o +call_mode.o : $(MODES) +$(MODES) io.o: atoms.o +$(MODES) main.o : io.o +testfuncs.o elements.o mode_create.o: subroutines.o +testfuncs.o : build_subroutines.o diff --git a/src/atoms.f90 b/src/atoms.f90 new file mode 100644 index 0000000..6da9189 --- /dev/null +++ b/src/atoms.f90 @@ -0,0 +1,1167 @@ +MODULE atoms +! +!********************************************************************************** +!* ATOMS * +!********************************************************************************** +!* This module contains subroutines concerning atoms. * +!********************************************************************************** +!* (C) Feb. 2014 - Pierre Hirel * +!* Université de Lille, Sciences et Technologies * +!* UMR CNRS 8207, UMET - C6, F-59655 Villeneuve D'Ascq, France * +!* pierre.hirel@univ-lille.fr * +!* Last modification: P. Hirel - 13 March 2018 * +!********************************************************************************** +!* This program is free software: you can redistribute it and/or modify * +!* it under the terms of the GNU General Public License as published by * +!* the Free Software Foundation, either version 3 of the License, or * +!* (at your option) any later version. * +!* * +!* This program is distributed in the hope that it will be useful, * +!* but WITHOUT ANY WARRANTY; without even the implied warranty of * +!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +!* GNU General Public License for more details. * +!* * +!* You should have received a copy of the GNU General Public License * +!* along with this program. If not, see . * +!********************************************************************************** +!* List of subroutines in this module: * +!* ATOMNUMBER gives the atomic number of an atom provided its symbol * +!* ATOMMASS gives the mass of an atom provided its symbol * +!* ATOMSPECIES gives an atom symbol provided its atomic number * +!* ATOMMASSSPECIES gives an atom symbol provided its mass * +!********************************************************************************** +! +! +! +use parameters +use functions + +INTEGER,PARAMETER,PUBLIC :: ATOMMAXZ = 118 !maximum value of atomic number +! +! +CONTAINS +! +!******************************************************** +! ATOMNUMBER +! This subroutine sets the atomic number of an atom +! depending on its species. +!******************************************************** +! +SUBROUTINE ATOMNUMBER(species,snumber) +! +IMPLICIT NONE +CHARACTER(LEN=2):: species, species2 +REAL(dp),INTENT(OUT):: snumber +! +! +!Make sure that the first letter is uppercase, the second one lowercase +species2(1:1) = StrUpCase(species(1:1)) +species2(2:2) = StrDnCase(species(2:2)) +! +SELECT CASE(species2) +! n=1 +CASE('H','D') + snumber=1.d0 +CASE('He') + snumber=2.d0 +! +! n=2 +CASE('Li') + snumber=3.d0 +CASE('Be') + snumber=4.d0 +CASE('B') + snumber=5.d0 +CASE('C') + snumber=6.d0 +CASE('N') + snumber=7.d0 +CASE('O') + snumber=8.d0 +CASE('F') + snumber=9.d0 +CASE('Ne') + snumber=10.d0 +! +! n=3 +CASE('Na') + snumber=11.d0 +CASE('Mg') + snumber=12.d0 +CASE('Al') + snumber=13.d0 +CASE('Si') + snumber=14.d0 +CASE('P') + snumber=15.d0 +CASE('S') + snumber=16.d0 +CASE('Cl') + snumber=17.d0 +CASE('Ar') + snumber=18.d0 +! +! n=4 +CASE('K') + snumber=19.d0 +CASE('Ca') + snumber=20.d0 +CASE('Sc') + snumber=21.d0 +CASE('Ti') + snumber=22.d0 +CASE('V') + snumber=23.d0 +CASE('Cr') + snumber=24.d0 +CASE('Mn') + snumber=25.d0 +CASE('Fe') + snumber=26.d0 +CASE('Co') + snumber=27.d0 +CASE('Ni') + snumber=28.d0 +CASE('Cu') + snumber=29.d0 +CASE('Zn') + snumber=30.d0 +CASE('Ga') + snumber=31.d0 +CASE('Ge') + snumber=32.d0 +CASE('As') + snumber=33.d0 +CASE('Se') + snumber=34.d0 +CASE('Br') + snumber=35.d0 +CASE('Kr') + snumber=36.d0 +! +! n=5 +CASE('Rb') + snumber=37.d0 +CASE('Sr') + snumber=38.d0 +CASE('Y') + snumber=39.d0 +CASE('Zr') + snumber=40.d0 +CASE('Nb') + snumber=41.d0 +CASE('Mo') + snumber=42.d0 +CASE('Tc') + snumber=43.d0 +CASE('Ru') + snumber=44.d0 +CASE('Rh') + snumber=45.d0 +CASE('Pd') + snumber=46.d0 +CASE('Ag') + snumber=47.d0 +CASE('Cd') + snumber=48.d0 +CASE('In') + snumber=49.d0 +CASE('Sn') + snumber=50.d0 +CASE('Sb') + snumber=51.d0 +CASE('Te') + snumber=52.d0 +CASE('I') + snumber=53.d0 +CASE('Xe') + snumber=54.d0 +! +! n=6 +CASE('Cs') + snumber=55.d0 +CASE('Ba') + snumber=56.d0 +! Lanthanides +CASE('La') + snumber=57.d0 +CASE('Ce') + snumber=58.d0 +CASE('Pr') + snumber=59.d0 +CASE('Nd') + snumber=60.d0 +CASE('Pm') + snumber=61.d0 +CASE('Sm') + snumber=62.d0 +CASE('Eu') + snumber=63.d0 +CASE('Gd') + snumber=64.d0 +CASE('Tb') + snumber=65.d0 +CASE('Dy') + snumber=66.d0 +CASE('Ho') + snumber=67.d0 +CASE('Er') + snumber=68.d0 +CASE('Tm') + snumber=69.d0 +CASE('Yb') + snumber=70.d0 +CASE('Lu') + snumber=71.d0 +! End of Lanthanides +CASE('Hf') + snumber=72.d0 +CASE('Ta') + snumber=73.d0 +CASE('W') + snumber=74.d0 +CASE('Re') + snumber=75.d0 +CASE('Os') + snumber=76.d0 +CASE('Ir') + snumber=77.d0 +CASE('Pt') + snumber=78.d0 +CASE('Au') + snumber=79.d0 +CASE('Hg') + snumber=80.d0 +CASE('Tl') + snumber=81.d0 +CASE('Pb') + snumber=82.d0 +CASE('Bi') + snumber=83.d0 +CASE('Po') + snumber=84.d0 +CASE('At') + snumber=85.d0 +CASE('Rn') + snumber=86.d0 +! +! n=7 +CASE('Fr') + snumber=87.d0 +CASE('Ra') + snumber=88.d0 +! Actinides +CASE('Ac') + snumber=89.d0 +CASE('Th') + snumber=90.d0 +CASE('Pa') + snumber=91.d0 +CASE('U') + snumber=92.d0 +CASE('Np') + snumber=93.d0 +CASE('Pu') + snumber=94.d0 +CASE('Am') + snumber=95.d0 +CASE('Cm') + snumber=96.d0 +CASE('Bk') + snumber=97.d0 +CASE('Cf') + snumber=98.d0 +CASE('Es') + snumber=99.d0 +CASE('Fm') + snumber=100.d0 +CASE('Md') + snumber=101.d0 +CASE('No') + snumber=102.d0 +CASE('Lr') + snumber=103.d0 +! End of actinides +CASE('Rf') + snumber=104.d0 +CASE('Db') + snumber=105.d0 +CASE('Sg') + snumber=106.d0 +CASE('Bh') + snumber=107.d0 +CASE('Hs') + snumber=108.d0 +CASE('Mt') + snumber=109.d0 +CASE('Ds') + snumber=110.d0 +CASE('Rg') + snumber=111.d0 +CASE('Cn') + snumber=112.d0 +CASE('Nh') + snumber=113.d0 +CASE('Fl') + snumber=114.d0 +CASE('Mc') + snumber=115.d0 +CASE('Lv') + snumber=116.d0 +CASE('Ts') + snumber=117.d0 +CASE('Og') + snumber=118.d0 +! +CASE DEFAULT + !If the species is not recognized + snumber=0.d0 +END SELECT +! +END SUBROUTINE ATOMNUMBER +! +! +! +!******************************************************** +! ATOMMASS +! This subroutine sets the mass of an atom +! depending on its species. The masses are from the +! National Institute of Standards and Technology (NIST): +! http://www.nist.gov/pml/data/periodic.cfm +!******************************************************** +! +SUBROUTINE ATOMMASS(species,smass) +! +IMPLICIT NONE +CHARACTER(LEN=2),INTENT(IN):: species +CHARACTER(LEN=2):: species2 +REAL(dp):: smass +! +! +!Make sure that the first letter is uppercase, the second one lowercase +species2(1:1) = StrUpCase(species(1:1)) +species2(2:2) = StrDnCase(species(2:2)) +! +! +SELECT CASE(species2) +! n=1 +CASE('H') + smass=1.008d0 +CASE('D') + smass=2.014101777d0 +CASE('He') + smass=4.002602d0 +! +! n=2 +CASE('Li') + smass=6.94d0 +CASE('Be') + smass=9.012182d0 +CASE('B') + smass=10.81d0 +CASE('C') + smass=12.011d0 +CASE('N') + smass=14.007d0 +CASE('O') + smass=15.999d0 +CASE('F') + smass=18.9984032d0 +CASE('Ne') + smass=20.1797d0 +! +! n=3 +CASE('Na') + smass=22.98976928d0 +CASE('Mg') + smass=24.305d0 +CASE('Al') + smass=26.9815386d0 +CASE('Si') + smass=28.085d0 +CASE('P') + smass=30.973762d0 +CASE('S') + smass=32.06d0 +CASE('Cl') + smass=35.45d0 +CASE('Ar') + smass=39.948d0 +! +! n=4 +CASE('K') + smass=39.0983d0 +CASE('Ca') + smass=40.078d0 +CASE('Sc') + smass=44.955912d0 +CASE('Ti') + smass=47.867d0 +CASE('V') + smass=50.9415d0 +CASE('Cr') + smass=51.9961d0 +CASE('Mn') + smass=54.938045d0 +CASE('Fe') + smass=55.845d0 +CASE('Co') + smass=58.933195d0 +CASE('Ni') + smass=58.6934d0 +CASE('Cu') + smass=63.546d0 +CASE('Zn') + smass=65.38d0 +CASE('Ga') + smass=69.723d0 +CASE('Ge') + smass=72.63d0 +CASE('As') + smass=74.9216d0 +CASE('Se') + smass=78.96d0 +CASE('Br') + smass=79.904d0 +CASE('Kr') + smass=83.798d0 +! +! n=5 +CASE('Rb') + smass=85.4678d0 +CASE('Sr') + smass=87.62d0 +CASE('Y') + smass=88.90585d0 +CASE('Zr') + smass=91.224d0 +CASE('Nb') + smass=92.90638d0 +CASE('Mo') + smass=95.96d0 +CASE('Tc') + smass=98.906d0 +CASE('Ru') + smass=101.07d0 +CASE('Rh') + smass=102.90550d0 +CASE('Pd') + smass=106.42d0 +CASE('Ag') + smass=107.8682d0 +CASE('Cd') + smass=112.411d0 +CASE('In') + smass=114.818d0 +CASE('Sn') + smass=118.71d0 +CASE('Sb') + smass=121.76d0 +CASE('Te') + smass=127.60d0 +CASE('I') + smass=126.90447d0 +CASE('Xe') + smass=131.293d0 +! +! n=6 +CASE('Cs') + smass=132.9054519d0 +CASE('Ba') + smass=137.327d0 +! Lanthanides +CASE('La') + smass=138.90547d0 +CASE('Ce') + smass=140.116d0 +CASE('Pr') + smass=140.90765d0 +CASE('Nd') + smass=144.242d0 +CASE('Pm') + smass=144.91d0 +CASE('Sm') + smass=150.36d0 +CASE('Eu') + smass=151.964d0 +CASE('Gd') + smass=157.25d0 +CASE('Tb') + smass=158.92535d0 +CASE('Dy') + smass=162.50d0 +CASE('Ho') + smass=164.93032d0 +CASE('Er') + smass=167.259d0 +CASE('Tm') + smass=168.93421d0 +CASE('Yb') + smass=173.054d0 +CASE('Lu') + smass=174.9668d0 +! End of Lanthanides +CASE('Hf') + smass=178.49d0 +CASE('Ta') + smass=180.94788d0 +CASE('W') + smass=183.84d0 +CASE('Re') + smass=186.207d0 +CASE('Os') + smass=190.23d0 +CASE('Ir') + smass=192.217d0 +CASE('Pt') + smass=195.084d0 +CASE('Au') + smass=196.966569d0 +CASE('Hg') + smass=200.59d0 +CASE('Tl') + smass=204.38d0 +CASE('Pb') + smass=207.2d0 +CASE('Bi') + smass=208.9804d0 +CASE('Po') + smass=209.98d0 +CASE('At') + smass=209.99d0 +CASE('Rn') + smass=222.02d0 +! +! n=7 +CASE('Fr') + smass=233.d0 +CASE('Ra') + smass=226.d0 +! Actinides +CASE('Ac') + smass=227.d0 +CASE('Th') + smass=232.03806d0 +CASE('Pa') + smass=231.03588d0 +CASE('U') + smass=238.02891d0 +CASE('Np') + smass=237.d0 +CASE('Pu') + smass=244.d0 +CASE('Am') + smass=243.d0 +CASE('Cm') + smass=247.d0 +CASE('Bk') + smass=247.d0 +CASE('Cf') + smass=251.d0 +CASE('Es') + smass=252.d0 +CASE('Fm') + smass=257.d0 +CASE('Md') + smass=258.d0 +CASE('No') + smass=259.d0 +CASE('Lr') + smass=262.d0 +! End of actinides +CASE('Rf') + smass=265.d0 +CASE('Db') + smass=268.d0 +CASE('Sg') + smass=271.d0 +CASE('Bh') + smass=270.d0 +CASE('Hs') + smass=277.d0 +CASE('Mt') + smass=276.d0 +CASE('Ds') + smass=281.d0 +CASE('Rg') + smass=280.d0 +CASE('Cn') + smass=285.17d0 +CASE('Uu') + smass=284.d0 +CASE('Fl') + smass=289.d0 +CASE('Mc') + smass=288.d0 +CASE('Lv') + smass=293.d0 +CASE('Ts') + smass=294.d0 +CASE('Og') + smass=294.d0 +! +CASE DEFAULT + !If the species is not recognized + smass=0.d0 +END SELECT +! +END SUBROUTINE ATOMMASS +! +! +! +!******************************************************** +! ATOMSPECIES +! This subroutine sets the species of an atom +! depending on its atomic number +!******************************************************** +! +SUBROUTINE ATOMSPECIES(snumber,species) +! +IMPLICIT NONE +REAL(dp),INTENT(IN):: snumber +CHARACTER(LEN=2):: species +INTEGER:: smint +! +smint = INT(snumber) +! +SELECT CASE(smint) +! n=1 +CASE(1) + species='H' +CASE(2) + species='He' +! +! n=2 +CASE(3) + species='Li' +CASE(4) + species='Be' +CASE(5) + species='B' +CASE(6) + species='C' +CASE(7) + species='N' +CASE(8) + species='O' +CASE(9) + species='F' +CASE(10) + species='Ne' +! +! n=3 +CASE(11) + species='Na' +CASE(12) + species='Mg' +CASE(13) + species='Al' +CASE(14) + species='Si' +CASE(15) + species='P' +CASE(16) + species='S' +CASE(17) + species='Cl' +CASE(18) + species='Ar' +! +! n=4 +CASE(19) + species='K' +CASE(20) + species='Ca' +CASE(21) + species='Sc' +CASE(22) + species='Ti' +CASE(23) + species='V' +CASE(24) + species='Cr' +CASE(25) + species='Mn' +CASE(26) + species='Fe' +CASE(27) + species='Co' +CASE(28) + species='Ni' +CASE(29) + species='Cu' +CASE(30) + species='Zn' +CASE(31) + species='Ga' +CASE(32) + species='Ge' +CASE(33) + species='As' +CASE(34) + species='Se' +CASE(35) + species='Br' +CASE(36) + species='Kr' +! +! n=5 +CASE(37) + species='Rb' +CASE(38) + species='Sr' +CASE(39) + species='Y' +CASE(40) + species='Zr' +CASE(41) + species='Nb' +CASE(42) + species='Mo' +CASE(43) + species='Tc' +CASE(44) + species='Ru' +CASE(45) + species='Rh' +CASE(46) + species='Pd' +CASE(47) + species='Ag' +CASE(48) + species='Cd' +CASE(49) + species='In' +CASE(50) + species='Sn' +CASE(51) + species='Sb' +CASE(52) + species='Te' +CASE(53) + species='I' +CASE(54) + species='Xe' +! +! n=6 +CASE(55) + species='Cs' +CASE(56) + species='Ba' +! Lanthanides +CASE(57) + species='La' +CASE(58) + species='Ce' +CASE(59) + species='Pr' +CASE(60) + species='Nd' +CASE(61) + species='Pm' +CASE(62) + species='Sm' +CASE(63) + species='Eu' +CASE(64) + species='Gd' +CASE(65) + species='Tb' +CASE(66) + species='Dy' +CASE(67) + species='Ho' +CASE(68) + species='Er' +CASE(69) + species='Tm' +CASE(70) + species='Yb' +CASE(71) + species='Lu' +! End of lanthanides +CASE(72) + species='Hf' +CASE(73) + species='Ta' +CASE(74) + species='W' +CASE(75) + species='Re' +CASE(76) + species='Os' +CASE(77) + species='Ir' +CASE(78) + species='Pt' +CASE(79) + species='Au' +CASE(80) + species='Hg' +CASE(81) + species='Tl' +CASE(82) + species='Pb' +CASE(83) + species='Bi' +CASE(84) + species='Po' +CASE(85) + species='At' +CASE(86) + species='Rn' +! +! n=7 +CASE(87) + species='Fr' +CASE(88) + species='Ra' +! Actinides +CASE(89) + species='Ac' +CASE(90) + species='Th' +CASE(91) + species='Pa' +CASE(92) + species='U' +CASE(93) + species='Np' +CASE(94) + species='Pu' +CASE(95) + species='Am' +CASE(96) + species='Cm' +CASE(97) + species='Bk' +CASE(98) + species='Cf' +CASE(99) + species='Es' +CASE(100) + species='Fm' +CASE(101) + species='Md' +CASE(102) + species='No' +CASE(103) + species='Lr' +! End of actinides +CASE(104) + species='Rf' +CASE(105) + species='Db' +CASE(106) + species='Sg' +CASE(107) + species='Bh' +CASE(108) + species='Hs' +CASE(109) + species='Mt' +CASE(110) + species='Ds' +CASE(111) + species='Rg' +CASE(112) + species='Cn' +CASE(113) + species='Nh' +CASE(114) + species='Fl' +CASE(115) + species='Mc' +CASE(116) + species='Lv' +CASE(117) + species='Ts' +CASE(118) + species='Og' +! +CASE DEFAULT + species='XX' +END SELECT +! +END SUBROUTINE ATOMSPECIES +! +! +! +!******************************************************** +! ATOMMASSSPECIES +! This subroutine sets the species of an atom, +! depending on its mass. The provided mass is rounded +! to the nearest integer to avoid numerical imprecision. +!******************************************************** +! +SUBROUTINE ATOMMASSSPECIES(smass,species) +! +IMPLICIT NONE +CHARACTER(LEN=2),INTENT(OUT):: species +INTEGER:: mass +REAL(dp),INTENT(IN):: smass +! +! +!Round mass to nearest integer +mass = NINT(smass) +! +! +SELECT CASE(mass) +! n=1 +CASE(1) + species = 'H' +CASE(2) + species = 'D' +CASE(4) + species = 'He' +! +! n=2 +CASE(6,7) + species = 'Li' +CASE(9) + species = 'Be' +CASE(10,11) + species = 'B' +CASE(12) + species = 'C' +CASE(14) + species = 'N' +CASE(16) + species = 'O' +CASE(19) + species = 'F' +CASE(20) + species = 'Ne' +! +! n=3 +CASE(22,23) + species = 'Na' +CASE(24,25) + species = 'Mg' +CASE(26,27) + species = 'Al' +CASE(28) + species = 'Si' +CASE(30,31) + species = 'P' +CASE(32) + species = 'S' +CASE(35,36) + species = 'Cl' +! +! n=4 +CASE(39) + species = 'K' +CASE(40) + species = 'Ca' !could also be Ar +CASE(44,45) + species = 'Sc' +CASE(47,48) + species = 'Ti' +CASE(50,51) + species = 'V' +CASE(52) + species = 'Cr' +CASE(54,55) + species = 'Mn' +CASE(56) + species = 'Fe' +CASE(58) + species = 'Ni' +CASE(59) + species = 'Co' +CASE(63,64) + species = 'Cu' +CASE(65,66) + species = 'Zn' +CASE(69,70) + species = 'Ga' +CASE(72,73) + species = 'Ge' +CASE(74,75) + species = 'As' +CASE(78,79) + species = 'Se' +CASE(80) + species = 'Br' +CASE(83,84) + species = 'Kr' +! +! n=5 +CASE(85,86) + species = 'Rb' +CASE(87,88) + species = 'Sr' +CASE(89) + species = 'Y' +CASE(91) + species = 'Zr' +CASE(93) + species = 'Nb' +CASE(95,96) + species = 'Mo' +CASE(98,99) + species = 'Tc' +CASE(101) + species = 'Ru' +CASE(102,103) + species = 'Rh' +CASE(106,107) + species = 'Pd' +CASE(108) + species = 'Ag' +CASE(112,113) + species = 'Cd' +CASE(114,115) + species = 'In' +CASE(118,119) + species = 'Sn' +CASE(121,122) + species = 'Sb' +CASE(128) + species = 'Te' +CASE(126,127) + species = 'I' +CASE(131) + species = 'Xe' +! +! n=6 +CASE(133) + species = 'Cs' +CASE(137) + species = 'Ba' +! Lanthanides +CASE(138,139) + species = 'La' +CASE(140) + species = 'Ce' +CASE(141) + species = 'Pr' +CASE(144) + species = 'Nd' +CASE(145) + species = 'Pm' +CASE(150) + species = 'Sm' +CASE(151,152) + species = 'Eu' +CASE(157) + species = 'Gd' +CASE(158,159) + species = 'Tb' +CASE(162,163) + species = 'Dy' +CASE(164,165) + species = 'Ho' +CASE(167) + species = 'Er' +CASE(168,169) + species = 'Tm' +CASE(173) + species = 'Yb' +CASE(175) + species = 'Lu' +! End of Lanthanides +CASE(178,179) + species = 'Hf' +CASE(180,181) + species = 'Ta' +CASE(183,184) + species = 'W' +CASE(186) + species = 'Re' +CASE(190) + species = 'Os' +CASE(192) + species = 'Ir' +CASE(195) + species = 'Pt' +CASE(196,197) + species = 'Au' +CASE(200,201) + species = 'Hg' +CASE(204) + species = 'Tl' +CASE(207) + species = 'Pb' +CASE(208,209) + species = 'Bi' +CASE(210) + species = 'Po' +CASE(222) + species = 'Rn' +! +! n=7 +CASE(233) + species = 'Fr' +CASE(226) + species = 'Ra' +! Actinides +CASE(227) + species = 'Ac' +CASE(232) + species = 'Th' +CASE(231) + species = 'Pa' +CASE(238) + species = 'U' +CASE(237) + species = 'Np' +CASE(244) + species = 'Pu' +CASE(243) + species = 'Am' +CASE(247) + species = 'Cm' +CASE(251) + species = 'Cf' +CASE(252) + species = 'Es' +CASE(257) + species = 'Fm' +CASE(258) + species = 'Md' +CASE(259) + species = 'No' +CASE(262) + species = 'Lr' +! End of actinides +CASE(265) + species = 'Rf' +CASE(268) + species = 'Db' +CASE(271) + species = 'Sg' +CASE(270) + species = 'Bh' +CASE(277) + species = 'Hs' +CASE(276) + species = 'Mt' +CASE(281) + species = 'Ds' +CASE(280) + species = 'Rg' +CASE(285) + species = 'Cn' +CASE(284) + species = 'Uu' +CASE(289) + species = 'Fl' +CASE(288) + species = 'Mc' +CASE(293) + species = 'Lv' +CASE(294) + species = 'Ts' +! +CASE DEFAULT + !If the species is not recognized + species = "Xx" +END SELECT +! +END SUBROUTINE ATOMMASSSPECIES +! +! +! +END MODULE atoms diff --git a/src/call_mode.f90 b/src/call_mode.f90 new file mode 100644 index 0000000..8fcdb0a --- /dev/null +++ b/src/call_mode.f90 @@ -0,0 +1,24 @@ +subroutine call_mode(arg_num,mode) + !This code is used to parse the command line argument for the mode information and calls the required + !mode module. + + use mode_create + use parameters + + implicit none + + integer, intent(in) :: arg_num + character(len=100), intent(in) :: mode + + select case(mode) + case('--create') + call create(arg_num) + + case default + print *, "Mode ", mode, " currently not accepted. Please check documentation for ", & + "accepted modes and rerun." + + stop 3 + + end select +end subroutine call_mode diff --git a/src/elements.f90 b/src/elements.f90 index 1fbacf6..53b827d 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -1,32 +1,187 @@ module elements + + !This module contains the elements data structures, structures needed for building regions + !and operations that are done on elements + use parameters + use subroutines + implicit none - use precision_comm_module + !Data structures used to represent the CAC elements. Each index represents an element + integer,allocatable :: tag_ele(:) !Element tag (used to keep track of id's + character(len=100), allocatable :: type_ele(:) !Element type + integer, allocatable :: size_ele(:), lat_ele(:) !Element siz + real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array - implicit none + integer :: ele_num=0 !Number of elements + + !Data structure used to represent atoms + integer, allocatable :: tag_atom(:) !atom id + character(len=100), allocatable:: type_atom(:) !atom type + real(kind =dp),allocatable :: r_atom(:,:) !atom position + integer :: atom_num=0 !Number of atoms + + !Variables for creating elements based on primitive cells + real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3) + integer :: cubic_faces(4,6) + + !Below are lattice type arrays which provide information on the general form of the elements. + !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased. + + integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type + + !These variables contain information on the basis, for simplicities sake we limit + !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point. + !This can be easily increased with no change to efficiency + integer :: max_basisnum, basisnum(10)!Max basis atom number, number of basis atoms in each lattice type + character(len=2) :: basis_type(10,10) !Atom type of each basis + real(kind=dp) :: basis_pos(3,10,10) !Basis atom positions + + !Simulation cell parameters + real(kind=dp) :: box_bd(6) + + + public + contains + + subroutine lattice_init + !This subroutine just intializes variables needed for building the different finite + !element types + + !First initialize the cubic cell + cubic_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, & + 1.0_dp, 0.0_dp, 0.0_dp, & + 1.0_dp, 1.0_dp, 0.0_dp, & + 0.0_dp, 1.0_dp, 0.0_dp, & + 0.0_dp, 0.0_dp, 1.0_dp, & + 1.0_dp, 0.0_dp, 1.0_dp, & + 1.0_dp, 1.0_dp, 1.0_dp, & + 0.0_dp, 1.0_dp, 1.0_dp /), & + shape(fcc_cell)) + + !Now we create a list containing the list of vertices needed to describe the 6 cube faces + cubic_faces(:,1) = (/ 1, 4, 8, 5 /) + cubic_faces(:,2) = (/ 2, 3, 7, 6 /) + cubic_faces(:,3) = (/ 1, 2, 6, 5 /) + cubic_faces(:,4) = (/ 3, 4, 8, 7 /) + cubic_faces(:,5) = (/ 1, 2, 3, 4 /) + cubic_faces(:,6) = (/ 5, 6, 7, 8 /) + + !!Now initialize the fcc primitive cell + fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, & + 0.5_dp, 0.5_dp, 0.0_dp, & + 0.5_dp, 1.0_dp, 0.5_dp, & + 0.0_dp, 0.5_dp, 0.5_dp, & + 0.5_dp, 0.0_dp, 0.5_dp, & + 1.0_dp, 0.5_dp, 0.5_dp, & + 1.0_dp, 1.0_dp, 1.0_dp, & + 0.5_dp, 0.5_dp, 1.0_dp /), & + shape(fcc_cell)) + + fcc_mat = reshape((/ 0.5_dp, 0.5_dp, 0.0_dp, & + 0.0_dp, 0.5_dp, 0.5_dp, & + 0.5_dp, 0.0_dp, 0.5_dp /), & + shape(fcc_mat)) + call matrix_inverse(fcc_mat,3,fcc_inv) + + max_basisnum = 0 + basisnum(:) = 0 + basis_pos(:,:,:) = 0.0_dp + ng_node(:) = 0 + end subroutine lattice_init + + subroutine cell_init(lapa,esize,ele_type, orient_mat, cell_mat) + !This subroutine uses the user provided information to transform the finite element cell to the correct + !size, orientation, and dimensions using the esize, lattice parameter, element_type, and orientation provided + !by the user + + integer, intent(in) :: esize + real(kind=dp), intent(in) :: lapa, orient_mat(3,3) + character(len=100), intent(in) :: ele_type + + real(kind=dp), dimension(3,max_ng_node), intent(out) :: cell_mat + + select case(trim(ele_type)) + + case('fcc') + cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, fcc_cell)) + case default + print *, "Element type ", trim(ele_type), " currently not accepted" + stop + end select + + end subroutine cell_init + + subroutine alloc_ele_arrays(n,m) + !This subroutine used to provide initial allocation for the atom and element arrays + + integer, intent(in) :: n, m !n-size of element arrays, m-size of atom arrays + + integer :: allostat + + !Allocate element arrays + if(n > 0) then + allocate(tag_ele(n), type_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), & + stat=allostat) + if(allostat > 0) then + print *, "Error allocating element arrays in elements.f90 because of: ", allostat + stop + end if + end if + + if(m > 0) then + !Allocate atom arrays + allocate(tag_atom(m), type_atom(m), r_atom(3,m), stat=allostat) + if(allostat > 0) then + print *, "Error allocating atom arrays in elements.f90 because of: ", allostat + stop + end if + end if + end subroutine + + subroutine add_element(type, size, lat, r) + !Subroutine which adds an element to the element arrays + integer, intent(in) :: size, lat + character(len=100), intent(in) :: type + real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node) + + ele_num = ele_num + 1 + tag_ele(ele_num) = ele_num + type_ele(ele_num) = type + size_ele(ele_num) = size + lat_ele(ele_num) = lat + r_node(:,:,:,ele_num) = r(:,:,:) + + + end subroutine add_element + + subroutine add_atom(type, r) + !Subroutine which adds an atom to the atom arrays + character(len=2), intent(in) :: type + real(kind=dp), intent(in), dimension(3) :: r - !This is the data structure which is used to represent the CAC elements - type element - integer :: tag = 0 !Element tag (used to keep track of id's - integer :: type = 0 !Lattice type of the element - integer :: size = 0 !Element size - !Nodal position array below only works for wedge or fcc elements - real(kind=wp) :: r_node(3,8) - end type + atom_num = atom_num+1 + tag_atom(atom_num) = atom_num + type_atom(atom_num) = type + r_atom(:,atom_num) = r(:) - !Finite element array - type(element), allocatable :: element_array(:) + end subroutine add_atom - integer :: ele_num + subroutine def_ng_node(n, element_types) + !This subroutine defines the maximum node number among n element types + integer, intent(in) :: n !Number of element types + character(len=100), dimension(n) :: element_types !Array of element type strings - !Data structure used to represent atoms - type atom - integer :: tag = 0 - integer :: type = 0 - real(kind =wp) :: r - end type + integer :: i - type(atom), allocatable :: atoms(:) + max_ng_node = 0 - integer :: atom_num + do i=1, n + select case(trim(adjustl(element_types(i)))) + case('fcc') + ng_node(i) = 8 + end select + if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i) + end do + end subroutine def_ng_node end module elements \ No newline at end of file diff --git a/src/functions.f90 b/src/functions.f90 new file mode 100644 index 0000000..e7c7a05 --- /dev/null +++ b/src/functions.f90 @@ -0,0 +1,235 @@ +module functions + + use parameters + + implicit none + + public + contains + + + ! Functions below this comment are taken from the functions module of atomsk + !******************************************************** + ! STRUPCASE + ! This function reads a string of any length + ! and capitalizes all letters. + !******************************************************** + FUNCTION StrUpCase (input_string) RESULT (UC_string) + ! + IMPLICIT NONE + CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz' + CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + CHARACTER(*),INTENT(IN):: input_string + CHARACTER(LEN(Input_String)):: UC_string !Upper-Case string that is produced + INTEGER:: i, n + ! + IF(LEN(input_string)==0) RETURN + UC_string = input_string + ! Loop over string elements + DO i=1,LEN(UC_string) + !Find location of letter in lower case constant string + n = INDEX( lower_case, UC_string(i:i) ) + !If current substring is a lower case letter, make it upper case + IF(n>0) THEN + UC_string(i:i) = upper_case(n:n) + ENDIF + END DO + ! + END FUNCTION StrUpCase + + !******************************************************** + ! STRDNCASE + ! This function reads a string of any length + ! and transforms all letters to lower case. + !******************************************************** + FUNCTION StrDnCase (input_string) RESULT (lc_string) + ! + IMPLICIT NONE + CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz' + CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + CHARACTER(*),INTENT(IN):: input_string + CHARACTER(LEN(Input_String)):: lc_string !Lower-Case string that is produced + INTEGER:: i, n + ! + IF(LEN(input_string)==0) RETURN + lc_string = input_string + ! Loop over string elements + DO i=1,LEN(lc_string) + !Find location of letter in upper case constant string + n = INDEX( upper_case, lc_string(i:i) ) + !If current substring is an upper case letter, make it lower case + IF(n>0) THEN + lc_string(i:i) = lower_case(n:n) + ENDIF + END DO +! +END FUNCTION StrDnCase + + pure function matrix_normal(a, n) + + integer :: i + integer, intent(in) :: n + real(kind = dp), dimension(n) :: v + real(kind = dp), dimension(n, n),intent(in) :: a + real(kind = dp), dimension(n,n) :: matrix_normal + + matrix_normal(:, :) = a(:, :) + + do i = 1, n + + v(:) = a(i,:) + + matrix_normal(i, :) = v(:) / norm2(v) + + end do + + return + end function matrix_normal + + pure function cross_product(a, b) + !Function which finds the cross product of two vectors + + real(kind = dp), dimension(3), intent(in) :: a, b + real(kind = dp), dimension(3) :: cross_product + + cross_product(1) = a(2) * b(3) - a(3) * b(2) + cross_product(2) = a(3) * b(1) - a(1) * b(3) + cross_product(3) = a(1) * b(2) - a(2) * b(1) + + return + end function cross_product + + pure function identity_mat(n) + !Returns the nxn identity matrix + + integer :: i + integer, intent(in) :: n + real(kind = dp), dimension(n, n) :: identity_mat + + identity_mat(:, :) = 0.0_dp + do i = 1, n + identity_mat(i, i) = 1.0_dp + end do + + return + end function identity_mat + + pure function triple_product(a, b, c) + !triple product between three 3*1 vectors + + real(kind = dp) :: triple_product + real(kind = dp), dimension(3), intent(in) :: a, b, c + triple_product = dot_product(a, cross_product(b, c)) + + return + end function triple_product + + function in_bd_lat(v, box_faces, box_norms) + !This function returns whether the point is within the transformed box boundaries. The transformed + !box being the transformed simulation cell in the lattice basis + + !Input/output variables + real(kind=dp), dimension(3), intent(in) :: v !integer lattice position + real(kind=dp), dimension(3,6), intent(in) :: box_faces !Centroid of all the box faces + real(kind=dp), dimension(3,6), intent(in) :: box_norms !Box face normals + logical :: in_bd_lat + + !Other variables + integer :: i + real(kind=dp) :: pt_to_face(3) + + in_bd_lat = .true. + + !Check if point is in box bounds, this works by comparing the dot product of the face normal and the + !vector from the point to the face. If the dot product is greater than 0 then the point is behind the face + !if it is equal to zero then the point is on the face, if is less than 0 the the point is in front of the face. + do i = 1, 6 + pt_to_face(:) = box_faces(:, i) - v + if(dot_product(pt_to_face, box_norms(:,i)) <= 0) then + in_bd_lat = .false. + exit + end if + end do + + return + end function in_bd_lat + + function in_block_bd(v, box_bd) + !This function returns whether a point is within a block in 3d + + !Input/output + real(kind=dp), dimension(3), intent(in) :: v + real(kind=dp), dimension(6), intent(in) :: box_bd + logical :: in_block_bd + + !Other variables + integer :: i + + in_block_bd = .true. + + do i =1 ,3 + !Check upper bound + if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then + in_block_bd =.false. + exit + !Check lower bound + else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then + in_block_bd = .false. + exit + end if + end do + end function in_block_bd + + function lcm(a,b) + !This function returns the smallest least common multiple of two numbers + + real(kind=dp), intent(in) :: a, b + real(kind=dp) :: lcm + + integer :: aint, bint, gcd, remainder, placeholder + + !Cast the vector positions to ints. There will be some error associated with this calculation + aint = a*10**2 + bint = b*10**2 + + !Calculate greated common divisor + gcd = aint + placeholder = bint + do while(placeholder /= 0) + remainder = modulo(gcd, placeholder) + gcd = placeholder + placeholder=remainder + end do + lcm = real((aint*bint),dp)/(real(gcd,dp))* 10.0_dp**(-2.0_dp) + end function lcm + + function is_neighbor(rl, rk, r_in, r_out) + !This function checks to see if two atoms are within a shell with an inner radius r_in and outer radius + !r_out + real(kind=dp), intent(in) :: r_in, r_out + real(kind=dp), dimension(3), intent(in) :: rl, rk + logical :: is_neighbor + + !Internal variable + real(kind=dp) :: rlk + + rlk = norm2(rk - rl) + is_neighbor=.true. + if((rlk>r_out).or.(rlk < r_in)) is_neighbor = .false. + + return + end function is_neighbor + + function is_equal(A, B) + !This function checks if too numbers are equal within a tolerance + real(kind=dp), intent(in) :: A, B + logical :: is_equal + + if((A>(B - 10.0_dp**-10)).and.(A < (B+10.0_dp**-10))) then + is_equal = .true. + else + is_equal = .false. + end if + return + end function is_equal +end module functions diff --git a/src/io.f90 b/src/io.f90 new file mode 100644 index 0000000..1bfa649 --- /dev/null +++ b/src/io.f90 @@ -0,0 +1,177 @@ +module io + + use elements + use parameters + use atoms + + implicit none + + integer :: outfilenum = 0 + character(len=100) :: outfiles(10) + + public + contains + + subroutine get_out_file(filename) + + implicit none + + character(len=100), intent(in) :: filename + character(len=100) :: temp_outfile + character(len=1) :: overwrite + logical :: file_exists + + !If no filename is provided then this function is called with none and prompts user input + if (filename=='none') then + print *, "Please specify a filename or extension to output to:" + read(*,*) temp_outfile + else + temp_outfile = filename + end if + + !Infinite loop which only exists if user provides valid filetype + overwrite = 'r' + do while(.true.) + + !Check to see if file exists, if it does then ask user if they would like to overwrite the file + inquire(file=trim(temp_outfile), exist=file_exists) + if (file_exists) then + if (overwrite == 'r') print *, "File ", trim(temp_outfile), " already exists. Would you like to overwrite? (Y/N)" + read(*,*) overwrite + if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then + print *, "Please specify a new filename with extension:" + read(*,*) temp_outfile + else if((scan(overwrite, "y") > 0).or.(scan(overwrite, "Y") > 0)) then + continue + else + print *, "Please pick either y or n" + read(*,*) overwrite + end if + + end if + + if (scan(temp_outfile,'.',.true.) == 0) then + print *, "No extension included on filename, please type a full filename that includes an extension." + read(*,*) temp_outfile + cycle + end if + select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) + case('xyz') + outfilenum=outfilenum+1 + outfiles(outfilenum) = temp_outfile + exit + case('lmp') + outfilenum=outfilenum+1 + outfiles(outfilenum) = temp_outfile + exit + + case default + print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", & + "please input a filename with extension from following list: xyz." + read(*,*) temp_outfile + + end select + end do + + end subroutine get_out_file + + + subroutine write_out + !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine + + integer :: i + + do i = 1, outfilenum + !Pull out the extension of the file and call the correct write subroutine + select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:)))) + case('xyz') + call write_xyz(outfiles(i)) + case('lmp') + call write_lmp(outfiles(i)) + case default + print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & + " is not accepted for writing. Please select from: xyz and try again" + stop + + end select + end do + end subroutine write_out + + + + subroutine write_xyz(file) + !This is the simplest visualization subroutine, it writes out all nodal positions and atom positions to an xyz file + character(len=100), intent(in) :: file + + integer :: node_num, i, inod, ibasis + + open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + + !Calculate total node number + node_num=0 + do i = 1, ele_num + node_num = node_num + basisnum(lat_ele(i))*ng_node(lat_ele(i)) + end do + + !Write total number of atoms + elements + write(11, '(i16)') node_num+atom_num + + !Write comment line + write(11, '(a)') "#Node + atom file created using cacmb" + + !Write nodal positions + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) + end do + end do + end do + + !Write atom positions + do i = 1, atom_num + write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i) + end do + + !Finish writing + close(11) + end subroutine write_xyz + + subroutine write_lmp(file) + + integer :: write_num, i + character(len=100), intent(in) :: file + !This subroutine writes out a .lmp style dump file + + open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + + !Comment line + write(11, '(a)') '# lmp file made with cacmb' + write(11, '(a)') + !Calculate total atom number + write_num = atom_num + !Write total number of atoms + elements + write(11, '(i16, a)') write_num, ' atoms' + !Write number of atom types + write(11, '(i16, a)') 1, ' atom types' + + write(11,'(a)') ' ' + !Write box bd + write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi' + write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi' + write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi' + + !Masses + write(11, '(a)') 'Masses' + write(11, '(a)') ' ' + write(11, '(i16, f23.15)') 1, 63.546 + write(11, '(a)') ' ' + + !Write atom positions + write(11, '(a)') 'Atoms' + write(11, '(a)') ' ' + do i = 1, atom_num + write(11, '(2i16, 3f23.15)') i, 1, r_atom(:,i) + end do + end subroutine write_lmp +end module io \ No newline at end of file diff --git a/src/lattice.f90 b/src/lattice.f90 deleted file mode 100644 index 4783a2b..0000000 --- a/src/lattice.f90 +++ /dev/null @@ -1,95 +0,0 @@ -module lattice - - use precision_comm_module - use subroutines - - implicit none - - integer :: atom_types - - !Atom type variables - character(len=2), dimension(10) :: atom_names - real(kind=wp), dimension(10) :: atom_masses - - !Lattice_type variables - integer :: lat_num - character(len=10), dimension(10) :: lattice_id, lattice_type - real(kind=wp), dimension(10) :: lapa - integer(kind=wp), dimension(10) :: basis_atom_num - integer(kind=wp), dimension(10,10) :: basis_type - real(kind=wp), dimension(3,10,10) ::basis_pos - - !Unit Cell variables - real(kind = wp) :: fcc_cell(3,8), fcc_mat(3,3) - - public - contains - - subroutine lattice_init - - !Initialize needed variables - lat_num=0 - basis_atom_num(:) = 0 - - !Initialize finite element cells to be used - - !First initialize the primitive fcc cell - fcc_cell = reshape((/ 0.0_wp, 0.0_wp, 0.0_wp, & - 0.5_wp, 0.5_wp, 0.0_wp, & - 0.5_wp, 1.0_wp, 0.5_wp, & - 0.0_wp, 0.5_wp, 0.5_wp, & - 0.5_wp, 0.0_wp, 0.5_wp, & - 1.0_wp, 0.5_wp, 0.5_wp, & - 1.0_wp, 1.0_wp, 1.0_wp, & - 0.5_wp, 0.5_wp, 1.0_wp /), & - shape(fcc_cell)) - - fcc_mat = reshape((/ 0.5_wp, 0.5_wp, 0.0_wp, & - 0.5_wp, 0.5_wp, 0.5_wp, & - 0.5_wp, 0.0_wp, 0.5_wp /), & - shape(fcc_mat)) - end subroutine lattice_init - !This subroutine defines the atom type arrays - subroutine atom_type_parse(line) - - character(len=100), intent(in) :: line - character(len=100) :: errorloc - integer :: ia, error - character(len=20) :: label - - read(line, *, iostat=error) label, atom_types, (atom_names(ia), atom_masses(ia), ia=1, atom_types) - errorloc="lattice:22" - call read_error_check(error,errorloc) - - end subroutine atom_type_parse - - !This subroutine defines the lattice types and the unit cells for the lattice types - subroutine lattice_parse(line) - character(len=100), intent(in) :: line - integer :: ia, error - character(len=20) :: label, kw - character(len=100) :: errorloc - - lat_num = lat_num + 1 - read(line, *, iostat=error) label, lattice_id(lat_num), lattice_type(lat_num), lapa(lat_num), kw - errorloc="lattice:77" - call read_error_check(error, errorloc) - - select case(kw) - case("type") - read(line(scan(line, "type"):), *, iostat=error) label, basis_type(1,1) - errorloc="lattice:56" - call read_error_check(error,errorloc) - case("basis") - read(line(scan(line, "basis"):), *, iostat=error) label, basis_atom_num(lat_num), (basis_type(ia, lat_num) ,& - basis_pos(1:3,ia,lat_num), ia = 1, basis_atom_num(lat_num)) - errorloc="lattice:59" - call read_error_check(error,errorloc) - case default - print *, "Keyword ", kw, " is not accepted in the lattice command" - stop "Exit with error" - end select - - - end subroutine lattice_parse -end module lattice \ No newline at end of file diff --git a/src/main.f90 b/src/main.f90 index 323f996..7f76b6b 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,52 +1,39 @@ program main - - use precision_comm_module - use elements - use lattice - use region + !**************************** CACmb ******************************* + !* CAC model building toolkit * + ! ____________ * + ! / / * + ! / / * + ! /___________/ * + ! _|_ _|_ _|____________ * + ! / / * + ! / / * + ! /___________/ * + ! * + !******************************************************************* - integer :: iosline, iospara - logical :: flags(4) - character(len=100) :: line, command, errorloc + use parameters + use elements + use io + - iosline = 0 - iospara = 0 - flags(:) = .false. + integer :: arg_num + character(len=100) :: mode call lattice_init - !Main command loop - do while (iosline == 0) + ! Command line parsing + arg_num = command_argument_count() - read(*, '(a)', iostat=iosline) line - errorloc="read_input:line" - call read_error_check(iosline, errorloc) + !Determine if a mode is being used and what it is. The first argument has to be the mode + !if a mode is being used + call get_command_argument(1, mode) - !Check for comment character (#) - if ((scan(line, '#')/= 1).and.(line/='')) then - read(line, *, iostat = iospara) command - errorloc="read_input:command" - call read_error_check(iosline, errorloc) + mode = trim(adjustl(mode)) + if (mode(1:2) == '--') then + call call_mode(arg_num, mode) + end if - select case(command) - case('atom_types') - call atom_type_parse(line) - flags(1) = .true. - case('lattice') - if(flags(1).eqv..false.) then - print *, "Please define atom types before defining lattice types" - stop 3 - end if - call lattice_parse(line) - flags(2) =.true. - ! case('region') - ! call build_region(line) - ! case('write') - ! call write_parse(line) - case default - print *, "The command ", trim(command), " is not currently accepted",& - " please check input script and try again." - end select - end if - end do + !Finish by writing the files + call write_out end program main \ No newline at end of file diff --git a/src/mode_create.f90 b/src/mode_create.f90 new file mode 100644 index 0000000..e6bb0be --- /dev/null +++ b/src/mode_create.f90 @@ -0,0 +1,432 @@ +module mode_create + !This mode is intended for creating element/atom regions and writing them to specific files. + + use parameters + use atoms + use io + use subroutines + + implicit none + + character(len=100) :: name, element_type + real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), & + orient_inv(3,3), box_vert(3,8), ox_bd(6), maxbd(3), lattice_space(3) + integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_num, lat_atom_num, bd_in_lat(6) + logical :: dup_flag, dim_flag + + real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) + public + contains + + subroutine create(arg_num) + ! Main subroutine which controls execution + + integer, intent(in) :: arg_num + character(len=100) :: textholder + + integer :: i, ibasis, inod + real(kind=dp) :: r(3), periodvone(3), periodvtwo(3) + real(kind=dp), allocatable :: r_node(:,:,:) + + !Initialize default parameters + orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient)) + cell_mat(:,:)=0.0_dp + name ='' + element_type = '' + esize=0 + lattice_parameter=0.0_dp + duplicate(:) = 0 + box_len(:) = 0.0_dp + dup_flag = .false. + dim_flag = .false. + basisnum = 0 + lat_num = 0 + lat_atom_num = 0 + + !First we parse the command + call parse_command(arg_num) + + ! Before building do a check on the file + if (outfilenum == 0) then + textholder = 'none' + call get_out_file(textholder) + end if + + !Now we setup the unit element and call other init subroutines + call def_ng_node(1, element_type) + + allocate(r_node(3,max_basisnum,max_ng_node)) + + if(dup_flag) then + + !We initialize the cell with a lattice_parameter of 1 because we will add the lattice parameter later + call cell_init(1.0_dp, esize, element_type, orient, cell_mat) + + + do i = 1, 8 + box_vert(:,i) = duplicate(:)*lattice_space(:)*cubic_cell(:,i) + origin(:) + end do + call matrix_inverse(orient,3,orient_inv) + !Now get the rotated box vertex positions in lattice space. Should be integer units + box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1 + !Find the new maxlen + maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2) + do i = 1, 3 + box_bd(2*i) = maxval(box_vert(i,:)) - 0.1_dp*lattice_space(i) + box_bd(2*i-1) = origin(i) + end do + + !and then call the build function with the correct transformation matrix + select case(trim(adjustl(element_type))) + case('fcc') + ! periodvone(:) = matmul(orient, reshape((/ 1, 1, 0 /),(/ 3 /))) + ! periodvtwo(:) = matmul(orient, reshape((/ 1, 1, 2 /),(/ 3 /))) + + ! do i = 1, 3 + ! if (periodvone(i) < lim_zero) then + ! ! box_bd(2*i) =floor(box_bd(2*i)/periodvtwo(i))*periodvtwo(i) + ! box_bd(2*i) = box_bd(2*i) - 0.5*periodvtwo(i) + ! else if(periodvtwo(i) < lim_zero) then + ! ! box_bd(2*i) =floor(box_bd(2*i)/periodvone(i))*periodvone(i) + ! box_bd(2*i) = box_bd(2*i) - 0.5*periodvone(i) + ! else + ! ! box_bd(2*i) = floor(box_bd(2*i)/lcm(periodvone(i),periodvtwo(i)))*lcm(periodvone(i),periodvtwo(i)) + ! box_bd(2*i) = box_bd(2*i) - 0.5*lcm(periodvone(i),periodvtwo(i)) + + ! end if + ! end do + + call lattice_in_box(box_lat_vert, fcc_mat) + case default + print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", & + "element type" + stop 3 + end select + + !Now that it is multiply by the lattice parameter + box_bd = box_bd*lattice_parameter + else if(dim_flag) then + continue + else + + call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) + !If the user doesn't pass any build instructions than we just put the cell mat into the element_array + call alloc_ele_arrays(1,0) + + !Add the basis atoms to the unit cell + do inod = 1, max_ng_node + do ibasis = 1, basisnum(1) + r_node(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) + end do + end do + + call add_element(element_type, esize, 1, r_node) + end if + + !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays + if(dup_flag.or.dim_flag) then + !Allocate variables + call alloc_ele_arrays(lat_num, lat_atom_num*basisnum(1)) + if(lat_atom_num > 0) then + !Check for periodicity + do i = 1, lat_atom_num + do ibasis = 1, basisnum(1) + call add_atom(basis_type(ibasis,1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis,1)) + end do + end do + deallocate(r_atom_lat) + end if + end if + + end subroutine create + !This subroutine parses the command and pulls out information needed for mode_create + subroutine parse_command(arg_num) + + integer, intent(in) :: arg_num + + integer :: arg_pos, ori_pos, i, j, arglen, stat + character(len=100) :: textholder + character(len=8) :: orient_string + + + !Pull out all required positional arguments + call get_command_argument(2, name, arglen) + if(arglen==0) STOP "Name is missing in mode create" + + call get_command_argument(3,element_type, arglen) + if(arglen==0) STOP "Element_type is missing in mode create" + + call get_command_argument(4,textholder, arglen) + if(arglen==0) STOP "Lattice Parameter is missing in mode create" + read(textholder, *, iostat=stat) lattice_parameter + if(stat > 0) STOP "Error reading lattice parameter" + + call get_command_argument(5, textholder, arglen) + if(arglen==0) STOP "Esize missing in mode create" + read(textholder, *, iostat=stat) esize + if(stat > 0) STOP "Error reading esize" + + arg_pos = 6 + !Check for optional keywords + do while(.true.) + if(arg_pos > command_argument_count()) exit + !Pull out the next argument which should either be a keyword or an option + call get_command_argument(arg_pos, textholder) + textholder=adjustl(textholder) + arg_pos=arg_pos+1 + + !Choose what to based on what the option string is + select case(trim(textholder)) + + !If orient command is passed extract the orientation to numeric array format + case('orient') + do i = 1, 3 + call get_command_argument(arg_pos, orient_string, arglen) + if (arglen==0) STOP "Missing orientation in orient command of mode create" + arg_pos = arg_pos+1 + ori_pos=2 + do j = 1,3 + if (orient_string(ori_pos:ori_pos) == '-') then + ori_pos = ori_pos + 1 + read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) + if (stat>0) STOP "Error reading orient value" + orient(i,j) = -orient(i,j) + ori_pos = ori_pos + 1 + else + read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) + if(stat>0) STOP "Error reading orient value" + ori_pos=ori_pos + 1 + end if + end do + end do + + + + !If the duplicate command is passed then we extract the information on the new bounds. + case('duplicate') + dup_flag = .true. + do i = 1, 3 + call get_command_argument(arg_pos, textholder) + read(textholder, *) duplicate(i) + arg_pos = arg_pos + 1 + end do + + case('origin') + do i = 1, 3 + call get_command_argument(arg_pos, textholder) + read(textholder, *) origin(i) + arg_pos = arg_pos + 1 + end do + print *, origin + !If a filetype is passed then we add name.ext to the outfiles list + case('xyz') + textholder = trim(adjustl(name)) //'.xyz' + call get_out_file(textholder) + + case default + !Check to see if it is an option command, if so then mode_create must be finished + if(textholder(1:1) == '-') then + exit + + !Check to see if a filename was passed + elseif(scan(textholder,'.',.true.) > 0) then + call get_out_file(textholder) + end if + end select + end do + + !Calculate the lattice periodicity length in lattice units + do i = 1, 3 + lattice_space(i) = norm2(orient(i,:)) + end do + + !Check special periodicity relations + select case(trim(adjustl(element_type))) + case('fcc') + do i = 1,3 + print *, orient(i,:) + !Check if one of the directions is 110 + if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.& + (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(orient(i,1),0.0_dp))).or.& + (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(orient(i,2),0.0_dp)))) then + + print *, '110', i + lattice_space(i) = 0.5_dp * lattice_space(i) + + !Check if one direction is 112 + else if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(abs(orient(i,3)),2.0_dp*abs(orient(i,1))))).or.& + (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(abs(orient(i,1)),2.0_dp*abs(orient(i,2))))).or.& + (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))& + then + + print *, '112 ', i + lattice_space(i) = 0.5_dp * lattice_space(i) + + end if + end do + end select + !Now normalize the orientation matrix + orient = matrix_normal(orient,3) + + !If we haven't defined a basis then define the basis and add the default basis atom type and position + if (basisnum(1) == 0) then + max_basisnum = 1 + basisnum(1) = 1 + basis_type(1,1) = name !If basis command not defined then we use name as the atom_name + basis_pos(:,1,1) = 0.0_dp + end if + end subroutine + + subroutine lattice_in_box(box_in_lat, transform_matrix) + !This subroutine returns all the lattice points in the box in r_lat + + !Inputs + integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space + real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space + !Internal variables + integer :: i, j, bd_in_lat(6), ix, iy, iz, numlatpoints + real(kind=dp) :: box_face_center(3,6), face_normals(3,6), Cx(2), Cy, Cz, A(2), v(3) + real(kind=dp), allocatable :: resize_lat_array(:,:) + + !First find the bounding lattice points (min and max points for the box in each dimension) + numlatpoints = 1 + do i = 1, 3 + bd_in_lat(2*i-1) = minval(box_in_lat(i,:)) + bd_in_lat(2*i) = maxval(box_in_lat(i,:)) + numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1)) + end do + + !Allocate the correct lat variable + select case(esize) + !Atomistics + case(2) + allocate(r_atom_lat(3,numlatpoints)) + case default + continue + end select + !Calculate the box_face centroids and box face normals. This is used in the centroid code. + box_face_center(:,:) = 0.0_dp + face_normals = reshape((/ -1.0_dp, 0.0_dp, 0.0_dp, & + 1.0_dp, 0.0_dp, 0.0_dp, & + 0.0_dp, -1.0_dp, 0.0_dp, & + 0.0_dp, 1.0_dp, 0.0_dp, & + 0.0_dp, 0.0_dp, -1.0_dp, & + 0.0_dp, 0.0_dp, 1.0_dp /),& + shape(face_normals)) + !Face normals + select case(trim(adjustl(element_type))) + case('fcc') + do i = 1, 6 + !Box face normal + face_normals(:,i) = matmul(fcc_inv, matmul(orient_inv, face_normals(:,i))) + end do + end select + + !Face centroids + do i =1, 6 + + !Initialize variables + A(:) = 0.0_dp + Cx(:) = 0.0_dp + Cy = 0.0_dp + Cz = 0.0_dp + + !Calculate all terms using a projection onto the xy and xz planes and then using the 2d equation + !for centroid of a plane. This is why we calculate the x centroid twice. + do j = 1, 4 + ! A(1) = A(1) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & + ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) + + ! !Centroid in x from xy plane + ! Cx(1) = Cx(1) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* & + ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & + ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) + + ! !Centroid in y from xy plane + ! Cy = Cy + (box_in_lat(2,cubic_faces(j,i))+box_in_lat(2,cubic_faces(j+1,i)))* & + ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & + ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) + + ! A(2) = A(2) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & + ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) + + ! !Centroid in x from xz plane + ! Cx(2) = Cx(2) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* & + ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & + ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) + + ! !Centroid in z from xz plane + ! Cz = Cz + (box_in_lat(3,cubic_faces(j,i))+box_in_lat(3,cubic_faces(j+1,i)))* & + ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & + ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) + + ! print *, j, i, Cx, Cy, Cz, A + Cx(1) = Cx(1) + box_in_lat(1,cubic_faces(j,i))*0.25 + Cy = Cy + box_in_lat(2,cubic_faces(j,i))*0.25 + Cz = Cz + box_in_lat(3,cubic_faces(j,i))*0.25 + + end do + + ! Cx = Cx * 1/(6*A) + ! if(Cx(1) /= Cx(2)) then + ! call error_message(7) + ! end if + ! Cy = Cy* 1/(6*A(1)) + ! Cz = Cz*1/(6*A(2)) + + box_face_center(:,i) = (/ Cx(1), Cy, Cz /) + end do + + !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case + !and one for the regular case + print *, box_bd + if (esize==2) then + !atomistics + do iz = bd_in_lat(5)-5, bd_in_lat(6)+5 + do iy = bd_in_lat(3)-5, bd_in_lat(4)+5 + do ix = bd_in_lat(1)-5, bd_in_lat(2)+5 + v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) + + !Transform point back to real space for easier checking + v = matmul(orient, matmul(transform_matrix,v)) + !If within the boundaries + if(in_block_bd(v, box_bd)) then + lat_atom_num = lat_atom_num + 1 + r_atom_lat(:,lat_atom_num) = v + end if + end do + end do + end do + else + continue + end if + + end subroutine lattice_in_box + + + subroutine error_message(errorid) + + integer, intent(in) :: errorid + + select case(errorid) + case(1) + STOP "Name is missing." + case(2) + print *, "Element_type is missing" + case(3) + print *, "Lattice_parameter is missing" + case(4) + print *, "Lattice parameter is not in float form" + case(5) + print *, "Esize is missing" + case(6) + print *, "Esize is not in integer form" + case(7) + print *, "Cx(1) should equal Cx(2) in plane centroid finding algorithm. Please double check implementation" + end select + + STOP 3 + end subroutine error_message + + +end module mode_create \ No newline at end of file diff --git a/src/parameters.f90 b/src/parameters.f90 new file mode 100644 index 0000000..0443622 --- /dev/null +++ b/src/parameters.f90 @@ -0,0 +1,8 @@ +module parameters + + implicit none + + integer, parameter :: dp= selected_real_kind(15,307) + real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), & + lim_large = huge(1.0_dp) +end module parameters \ No newline at end of file diff --git a/src/precision_comm_module.f90 b/src/precision_comm_module.f90 deleted file mode 100644 index 7532233..0000000 --- a/src/precision_comm_module.f90 +++ /dev/null @@ -1,13 +0,0 @@ - module precision_comm_module - - implicit none - - integer, parameter :: & - dp = selected_real_kind(15, 307), & ! double real - qp = selected_real_kind(33, 4931), & ! quadrupole real - wp = dp - - integer, save :: & - mpi_wp - - end module precision_comm_module diff --git a/src/region.f90 b/src/region.f90 deleted file mode 100644 index f1a9d33..0000000 --- a/src/region.f90 +++ /dev/null @@ -1,14 +0,0 @@ -module region - use precision_comm_module - - implicit none - - public - contains - - subroutine build_region(line) - - character(len=100), intent(in) :: line - - end subroutine build_region -end module region \ No newline at end of file diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 71aa475..0682037 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -1,13 +1,12 @@ module subroutines - - use precision_comm_module - + use parameters + use functions implicit none + integer :: allostat, deallostat public contains - !This subroutine is just used to break the code and exit on an error subroutine read_error_check(para, loc) @@ -18,7 +17,132 @@ module subroutines print *, "Read error in ", trim(loc), " because of ", para stop "Exit with error" end if + end subroutine + subroutine matrix_inverse(a, n, a_inv) + + integer :: i, j, k, piv_loc + integer, intent(in) :: n + real(kind = dp) :: coeff, sum_l, sum_u + real(kind = dp), dimension(n) :: b, x, y, b_piv + real(kind = dp), dimension(n, n) :: l, u, p + real(kind = dp), dimension(n, n), intent(in) :: a + real(kind = dp), dimension(n, n), intent(out) :: a_inv + real(kind = dp), allocatable :: v(:), u_temp(:), l_temp(:), p_temp(:) + + l(:, :) = identity_mat(n) + u(:, :) = a(:, :) + p(:, :) = identity_mat(n) + !LU decomposition with partial pivoting + do j = 1, n-1 + + allocate(v(n-j+1), stat = allostat) + if(allostat /=0 ) then + print *, 'Fail to allocate v in matrix_inverse' + stop + end if + + v(:) = u(j:n, j) + if(maxval(abs(v)) < lim_zero) then + print *, 'Fail to inverse matrix', a + stop + end if + + piv_loc = maxloc(abs(v), 1) + deallocate(v, stat = deallostat) + + if(deallostat /=0 ) then + print *, 'Fail to deallocate v in matrix_inverse' + stop + end if + !partial pivoting + if(piv_loc /= 1) then + + allocate( u_temp(n-j+1), p_temp(n), stat = allostat) + if(allostat /=0 ) then + print *, 'Fail to allocate p_temp and/or u_temp in matrix_inverse' + stop + end if + + u_temp(:) = u(j, j:n) + u(j, j:n) = u(piv_loc+j-1, j:n) + u(piv_loc+j-1, j:n) = u_temp(:) + p_temp(:) = p(j, :) + p(j, :) = p(piv_loc+j-1, :) + p(piv_loc+j-1, :) = p_temp(:) + + deallocate( u_temp, p_temp, stat = deallostat) + if(deallostat /=0 ) then + print *, 'Fail to deallocate p_temp and/or u_temp in matrix_inverse' + stop + end if + + if(j > 1) then + + allocate( l_temp(j-1), stat = allostat) + if(allostat /= 0) then + print *, 'Fail to allocate l_temp in matrix_inverse' + stop + end if + + l_temp(:) = l(j, 1:j-1) + l(j, 1:j-1) = l(piv_loc+j-1, 1:j-1) + l(piv_loc+j-1, 1:j-1) = l_temp(:) + + deallocate( l_temp, stat = deallostat) + if(deallostat /=0 ) then + print *, 'Fail to deallocate l_temp in matrix_inverse' + stop + end if + end if + end if + + !LU decomposition + do i = j+1, n + coeff = u(i, j)/u(j, j) + l(i, j) = coeff + u(i, j:n) = u(i, j:n)-coeff*u(j, j:n) + end do + + end do + + a_inv(:, :) = 0.0_dp + do j = 1, n + b(:) = 0.0_dp + b(j) = 1.0_dp + b_piv(:) = matmul(p, b) + !Now we have LUx = b_piv + !the first step is to solve y from Ly = b_piv + !forward substitution + do i = 1, n + if(i == 1) then + y(i) = b_piv(i)/l(i, i) + else + sum_l = 0 + do k = 1, i-1 + sum_l = sum_l+l(i, k)*y(k) + end do + y(i) = (b_piv(i)-sum_l)/l(i, i) + end if + end do + !then we solve x from ux = y + !backward subsitution + do i = n, 1, -1 + if(i == n) then + x(i) = y(i)/u(i, i) + else + sum_u = 0 + do k = i+1, n + sum_u = sum_u+u(i, k)*x(k) + end do + x(i) = (y(i)-sum_u)/u(i, i) + end if + end do + !put x into j column of a_inv + a_inv(:, j) = x(:) + end do + return + end subroutine matrix_inverse end module subroutines \ No newline at end of file From ff3fc5e20ab625e742a917d79e1b023c1aba04aa Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 27 Nov 2019 15:10:28 -0500 Subject: [PATCH 3/8] First working version of model builder with several output file types and mode_create working --- src/Makefile | 13 +- src/call_mode.f90 | 2 +- src/elements.f90 | 130 ++++++++++++++++++- src/functions.f90 | 2 +- src/io.f90 | 118 ++++++++++++++++-- src/mode_create.f90 | 297 ++++++++++++++++++++++++++------------------ 6 files changed, 419 insertions(+), 143 deletions(-) diff --git a/src/Makefile b/src/Makefile index c83df06..58507b4 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,24 +1,24 @@ FC=ifort -#FFLAGS=-c -mcmodel=large -debug -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -FFLAGS=-c -mcmodel=large -Ofast +FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone +#FFLAGS=-c -mcmodel=large -Ofast MODES=mode_create.o -OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o build_subroutines.o $(MODES) +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o $(MODES) .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o cacmb: $(OBJECTS) - $(FC) $(OBJECTS) -o $@ + $(FC) $(FFLAGS) $(OBJECTS) -o $@ .f90.o: - $(FC) $(FFLAGS) $< + $(FC) $(FFLAGS) -c $< .PHONY: clean clean: $(RM) cacmb *.o testfuncs: testfuncs.o functions.o subroutines.o - $(FC) testfuncs.o functions.o build_subroutines.o subroutines.o elements.o -o $@ + $(FC) testfuncs.o functions.o subroutines.o elements.o -o $@ .PHONY: cleantest cleantest: @@ -31,4 +31,3 @@ call_mode.o : $(MODES) $(MODES) io.o: atoms.o $(MODES) main.o : io.o testfuncs.o elements.o mode_create.o: subroutines.o -testfuncs.o : build_subroutines.o diff --git a/src/call_mode.f90 b/src/call_mode.f90 index 8fcdb0a..1ad4f0f 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -12,7 +12,7 @@ subroutine call_mode(arg_num,mode) select case(mode) case('--create') - call create(arg_num) + call create() case default print *, "Mode ", mode, " currently not accepted. Please check documentation for ", & diff --git a/src/elements.f90 b/src/elements.f90 index 53b827d..21de971 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -13,13 +13,17 @@ module elements real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array integer :: ele_num=0 !Number of elements + integer :: node_num=0 !Total number of nodes !Data structure used to represent atoms - integer, allocatable :: tag_atom(:) !atom id - character(len=100), allocatable:: type_atom(:) !atom type + integer, allocatable :: tag_atom(:), type_atom(:)!atom id real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms + !Mapping atom type to provided name + character(len=2), dimension(10) :: type_to_name + integer :: atom_types = 0 + !Variables for creating elements based on primitive cells real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3) integer :: cubic_faces(4,6) @@ -28,12 +32,12 @@ module elements !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased. integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type + integer :: max_esize=0 !Maximum number of atoms per side of element !These variables contain information on the basis, for simplicities sake we limit !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point. !This can be easily increased with no change to efficiency - integer :: max_basisnum, basisnum(10)!Max basis atom number, number of basis atoms in each lattice type - character(len=2) :: basis_type(10,10) !Atom type of each basis + integer :: max_basisnum, basisnum(10), basis_type(10,10)!Max basis atom number, number of basis atoms in each lattice type real(kind=dp) :: basis_pos(3,10,10) !Basis atom positions !Simulation cell parameters @@ -150,13 +154,14 @@ module elements size_ele(ele_num) = size lat_ele(ele_num) = lat r_node(:,:,:,ele_num) = r(:,:,:) + node_num = node_num + ng_node(lat) end subroutine add_element subroutine add_atom(type, r) !Subroutine which adds an atom to the atom arrays - character(len=2), intent(in) :: type + integer, intent(in) :: type real(kind=dp), intent(in), dimension(3) :: r atom_num = atom_num+1 @@ -166,6 +171,32 @@ module elements end subroutine add_atom + subroutine add_atom_type(type, inttype) + !This subroutine adds a new atom type to the module list of atom types + character(len=2), intent(in) :: type + integer, intent(out) :: inttype + + integer :: i + logical :: exists + + exists = .false. + do i=1, 10 + if(type == type_to_name(i)) exists = .true. + inttype = i + end do + + if (exists.eqv..false.) then + atom_types = atom_types+1 + if(atom_types > 10) then + print *, "Defined atom types are greater than 10 which is currently not supported." + stop 3 + end if + type_to_name(atom_types) = type + inttype = atom_types + end if + return + end subroutine add_atom_type + subroutine def_ng_node(n, element_types) !This subroutine defines the maximum node number among n element types integer, intent(in) :: n !Number of element types @@ -184,4 +215,93 @@ module elements if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i) end do end subroutine def_ng_node + + subroutine set_max_esize + !This subroutine sets the maximum esize + max_esize=maxval(size_ele) + end subroutine + + subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp) + !This subroutine returns the interpolated atoms from the elements. + + !Arguments + character(len=100), intent(in) :: type !The type of element that it is + integer, intent(in) :: esize !The number of atoms per side + integer, intent(in) :: lat_type !The integer lattice type of the element + real(kind=dp), dimension(3,max_basisnum, max_ng_node), intent(in) :: r_in !Nodal positions + integer, dimension(max_basisnum*max_esize**3), intent(out) :: type_interp !Interpolated atomic positions + real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions + + !Internal variables + integer :: i, it, is, ir, ibasis, inod, ia, bnum, lat_type_temp + real(kind=dp), allocatable :: a_shape(:) + real(kind=dp) :: t, s, r + + !Initialize some variables + r_interp(:,:) = 0.0_dp + type_interp(:) = 0 + ia = 0 + + !Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means + !the basis is 0,0,0, and the type doesn't matter + + select case(lat_type) + case(0) + bnum=1 + lat_type_temp = 1 + case default + bnum = basisnum(lat_type) + lat_type_temp = lat_type + end select + + select case(trim(adjustl(type))) + case('fcc') + allocate(a_shape(8)) + !Now loop over all the possible sites + do it = 1, esize + t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2) + do is =1, esize + s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2) + do ir = 1, esize + r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2) + call rhombshape(r,s,t,a_shape) + + do ibasis = 1, bnum + ia = ia + 1 + do inod = 1, 8 + type_interp(ia) = basis_type(ibasis,lat_type_temp) + r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod) + + end do + end do + + end do + end do + end do + end select + + if (ia /= esize**3) then + print *, "Incorrect interpolation" + stop 3 + end if + return + end subroutine interpolate_atoms + + subroutine rhombshape(r,s,t, shape_fun) + !Shape function for rhombohedral elements + real(kind=8), intent(in) :: r, s, t + real(kind=8), intent(out) :: shape_fun(8) + + shape_fun(1) = (1.0-r)*(1.0-s)*(1.0-t)/8.0 + shape_fun(2) = (1.0+r)*(1.0-s)*(1.0-t)/8.0 + shape_fun(3) = (1.0+r)*(1.0+s)*(1.0-t)/8.0 + shape_fun(4) = (1.0-r)*(1.0+s)*(1.0-t)/8.0 + shape_fun(5) = (1.0-r)*(1.0-s)*(1.0+t)/8.0 + shape_fun(6) = (1.0+r)*(1.0-s)*(1.0+t)/8.0 + shape_fun(7) = (1.0+r)*(1.0+s)*(1.0+t)/8.0 + shape_fun(8) = (1.0-r)*(1.0+s)*(1.0+t)/8.0 + + + return + end subroutine rhombshape end module elements \ No newline at end of file diff --git a/src/functions.f90 b/src/functions.f90 index e7c7a05..eb8e813 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -225,7 +225,7 @@ END FUNCTION StrDnCase real(kind=dp), intent(in) :: A, B logical :: is_equal - if((A>(B - 10.0_dp**-10)).and.(A < (B+10.0_dp**-10))) then + if((A>(B - 10.0_dp**(-10))).and.(A < (B+10.0_dp**(-10)))) then is_equal = .true. else is_equal = .false. diff --git a/src/io.f90 b/src/io.f90 index 1bfa649..7850ee2 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -64,10 +64,13 @@ module io outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit - + case('vtk') + outfilenum=outfilenum+1 + outfiles(outfilenum)=temp_outfile + exit case default print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", & - "please input a filename with extension from following list: xyz." + "please input a filename with extension from following list: xyz, lmp, vtk." read(*,*) temp_outfile end select @@ -81,6 +84,9 @@ module io integer :: i + !Find max esize which will be needed later + call set_max_esize + do i = 1, outfilenum !Pull out the extension of the file and call the correct write subroutine select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:)))) @@ -88,6 +94,8 @@ module io call write_xyz(outfiles(i)) case('lmp') call write_lmp(outfiles(i)) + case('vtk') + call write_vtk(outfiles(i)) case default print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & " is not accepted for writing. Please select from: xyz and try again" @@ -138,10 +146,10 @@ module io end subroutine write_xyz subroutine write_lmp(file) - - integer :: write_num, i - character(len=100), intent(in) :: file !This subroutine writes out a .lmp style dump file + character(len=100), intent(in) :: file + integer :: write_num, i, iatom, type_interp(max_basisnum*max_esize**3) + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), mass open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') @@ -150,10 +158,13 @@ module io write(11, '(a)') !Calculate total atom number write_num = atom_num + do i = 1,ele_num + if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3 + end do !Write total number of atoms + elements write(11, '(i16, a)') write_num, ' atoms' !Write number of atom types - write(11, '(i16, a)') 1, ' atom types' + write(11, '(i16, a)') atom_types, ' atom types' write(11,'(a)') ' ' !Write box bd @@ -163,15 +174,106 @@ module io !Masses write(11, '(a)') 'Masses' + write(11, '(a)') ' ' - write(11, '(i16, f23.15)') 1, 63.546 + do i =1, atom_types + call atommass(type_to_name(i),mass) + write(11, '(i16, f23.15)') i, mass + end do write(11, '(a)') ' ' !Write atom positions write(11, '(a)') 'Atoms' write(11, '(a)') ' ' do i = 1, atom_num - write(11, '(2i16, 3f23.15)') i, 1, r_atom(:,i) + write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) + end do + + !Write refined element atomic positions + do i = 1, ele_num + call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp) + select case(trim(adjustl(type_ele(i)))) + case('fcc') + do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom) + end do + end select end do end subroutine write_lmp + + subroutine write_vtk(file) + !This subroutine writes out a vtk style dump file + integer :: i, j, inod, ibasis + character(len=100), intent(in) :: file + +1 format('# vtk DataFile Version 4.0.1', / & + 'CAC output -- cg', / & + 'ASCII') +11 format('# vtk DataFile Version 4.0.1', / & + 'CACmb output -- atoms', / & + 'ASCII') +2 format('DATASET UNSTRUCTURED_GRID') +3 format('POINTS', i16, ' float') +4 format(/'CELLS', 2i16) +5 format(/'CELL_TYPES', i16) +12 format(/'CELL_DATA', i16) +16 format(/'POINT_DATA', i16) +17 format('SCALARS weight float', / & + 'LOOKUP_TABLE default') +18 format('SCALARS atom_type float', / & + 'LOOKUP_TABLE default') + +20 format('SCALARS lattice_type float', /& + 'LOOKUP_TABLE default') + + !First we write the vtk file containing the atoms + open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind') + + write(11, 11) + write(11, 2) + write(11, 3) atom_num + do i = 1, atom_num + write(11, '(3f23.15)') r_atom(:,i) + end do + write(11,4) atom_num, atom_num*2 + do i = 1, atom_num + write(11, '(2i16)') 1, i-1 + end do + write(11, 5) atom_num + do i = 1, atom_num + write(11, '(i16)') 1 + end do + write(11, 16) atom_num + write(11, 18) + do i = 1, atom_num + write(11, '(i16)') type_atom(i) + end do + close(11) + + open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind') + write(11,1) + write(11,2) + write(11,3) node_num + do i = 1, ele_num + do inod=1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + write(11, '(3f23.1)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i)) + end do + end do + end do + write(11, 4) ele_num, ele_num + node_num + do i =1, ele_num + write(11, '(9i16)') ng_node(lat_ele(i)), (j, j = (i-1)*ng_node(lat_ele(i)), i*ng_node(lat_ele(i))-1) + end do + write(11,5) ele_num + do i = 1, ele_num + if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12 + end do + write(11,12) ele_num + write(11,20) + do i = 1, ele_num + write(11, '(i16)') lat_ele(i) + end do + close(11) + end subroutine end module io \ No newline at end of file diff --git a/src/mode_create.f90 b/src/mode_create.f90 index e6bb0be..52342e2 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -5,28 +5,27 @@ module mode_create use atoms use io use subroutines + use elements implicit none character(len=100) :: name, element_type real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), & - orient_inv(3,3), box_vert(3,8), ox_bd(6), maxbd(3), lattice_space(3) - integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_num, lat_atom_num, bd_in_lat(6) + orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3) + integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6) logical :: dup_flag, dim_flag real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) public contains - subroutine create(arg_num) + subroutine create() ! Main subroutine which controls execution - integer, intent(in) :: arg_num character(len=100) :: textholder integer :: i, ibasis, inod - real(kind=dp) :: r(3), periodvone(3), periodvtwo(3) - real(kind=dp), allocatable :: r_node(:,:,:) + real(kind=dp), allocatable :: r_node_temp(:,:,:) !Initialize default parameters orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient)) @@ -40,11 +39,11 @@ module mode_create dup_flag = .false. dim_flag = .false. basisnum = 0 - lat_num = 0 + lat_ele_num = 0 lat_atom_num = 0 !First we parse the command - call parse_command(arg_num) + call parse_command() ! Before building do a check on the file if (outfilenum == 0) then @@ -55,7 +54,7 @@ module mode_create !Now we setup the unit element and call other init subroutines call def_ng_node(1, element_type) - allocate(r_node(3,max_basisnum,max_ng_node)) + allocate(r_node_temp(3,max_basisnum,max_ng_node)) if(dup_flag) then @@ -64,7 +63,7 @@ module mode_create do i = 1, 8 - box_vert(:,i) = duplicate(:)*lattice_space(:)*cubic_cell(:,i) + origin(:) + box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + origin(:) end do call matrix_inverse(orient,3,orient_inv) !Now get the rotated box vertex positions in lattice space. Should be integer units @@ -72,31 +71,14 @@ module mode_create !Find the new maxlen maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2) do i = 1, 3 - box_bd(2*i) = maxval(box_vert(i,:)) - 0.1_dp*lattice_space(i) - box_bd(2*i-1) = origin(i) + box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i) + box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i) end do - !and then call the build function with the correct transformation matrix select case(trim(adjustl(element_type))) case('fcc') - ! periodvone(:) = matmul(orient, reshape((/ 1, 1, 0 /),(/ 3 /))) - ! periodvtwo(:) = matmul(orient, reshape((/ 1, 1, 2 /),(/ 3 /))) - - ! do i = 1, 3 - ! if (periodvone(i) < lim_zero) then - ! ! box_bd(2*i) =floor(box_bd(2*i)/periodvtwo(i))*periodvtwo(i) - ! box_bd(2*i) = box_bd(2*i) - 0.5*periodvtwo(i) - ! else if(periodvtwo(i) < lim_zero) then - ! ! box_bd(2*i) =floor(box_bd(2*i)/periodvone(i))*periodvone(i) - ! box_bd(2*i) = box_bd(2*i) - 0.5*periodvone(i) - ! else - ! ! box_bd(2*i) = floor(box_bd(2*i)/lcm(periodvone(i),periodvtwo(i)))*lcm(periodvone(i),periodvtwo(i)) - ! box_bd(2*i) = box_bd(2*i) - 0.5*lcm(periodvone(i),periodvtwo(i)) - - ! end if - ! end do - - call lattice_in_box(box_lat_vert, fcc_mat) + + call build_with_rhomb(box_lat_vert, fcc_mat) case default print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", & "element type" @@ -116,33 +98,45 @@ module mode_create !Add the basis atoms to the unit cell do inod = 1, max_ng_node do ibasis = 1, basisnum(1) - r_node(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) + r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) end do end do - - call add_element(element_type, esize, 1, r_node) + do i = 1,3 + box_bd(2*i) = maxval(r_node_temp(i,:,:)) + box_bd(2*i-1) = origin(i) + end do + call add_element(element_type, esize, 1, r_node_temp) end if !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays if(dup_flag.or.dim_flag) then !Allocate variables - call alloc_ele_arrays(lat_num, lat_atom_num*basisnum(1)) + call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1)) if(lat_atom_num > 0) then - !Check for periodicity do i = 1, lat_atom_num - do ibasis = 1, basisnum(1) + do ibasis = 1, basisnum(1) call add_atom(basis_type(ibasis,1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis,1)) end do end do deallocate(r_atom_lat) end if + + if(lat_ele_num > 0) then + do i = 1, lat_ele_num + do inod= 1, ng_node(1) + do ibasis = 1, basisnum(1) + r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis,1) + end do + end do + call add_element(element_type, esize, 1, r_node_temp) + end do + end if end if end subroutine create !This subroutine parses the command and pulls out information needed for mode_create - subroutine parse_command(arg_num) + subroutine parse_command() - integer, intent(in) :: arg_num integer :: arg_pos, ori_pos, i, j, arglen, stat character(len=100) :: textholder @@ -217,7 +211,6 @@ module mode_create read(textholder, *) origin(i) arg_pos = arg_pos + 1 end do - print *, origin !If a filetype is passed then we add name.ext to the outfiles list case('xyz') textholder = trim(adjustl(name)) //'.xyz' @@ -244,13 +237,11 @@ module mode_create select case(trim(adjustl(element_type))) case('fcc') do i = 1,3 - print *, orient(i,:) !Check if one of the directions is 110 if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.& (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(orient(i,1),0.0_dp))).or.& (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(orient(i,2),0.0_dp)))) then - print *, '110', i lattice_space(i) = 0.5_dp * lattice_space(i) !Check if one direction is 112 @@ -259,7 +250,6 @@ module mode_create (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))& then - print *, '112 ', i lattice_space(i) = 0.5_dp * lattice_space(i) end if @@ -272,21 +262,27 @@ module mode_create if (basisnum(1) == 0) then max_basisnum = 1 basisnum(1) = 1 - basis_type(1,1) = name !If basis command not defined then we use name as the atom_name + call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name basis_pos(:,1,1) = 0.0_dp end if end subroutine - subroutine lattice_in_box(box_in_lat, transform_matrix) + subroutine build_with_rhomb(box_in_lat, transform_matrix) !This subroutine returns all the lattice points in the box in r_lat !Inputs integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space !Internal variables - integer :: i, j, bd_in_lat(6), ix, iy, iz, numlatpoints - real(kind=dp) :: box_face_center(3,6), face_normals(3,6), Cx(2), Cy, Cz, A(2), v(3) + integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, & + type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o + real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3) real(kind=dp), allocatable :: resize_lat_array(:,:) + logical, allocatable :: lat_points(:,:,:) + logical :: node_in_bd(8) + + !Do some value initialization + max_esize = esize !First find the bounding lattice points (min and max points for the box in each dimension) numlatpoints = 1 @@ -304,82 +300,10 @@ module mode_create case default continue end select - !Calculate the box_face centroids and box face normals. This is used in the centroid code. - box_face_center(:,:) = 0.0_dp - face_normals = reshape((/ -1.0_dp, 0.0_dp, 0.0_dp, & - 1.0_dp, 0.0_dp, 0.0_dp, & - 0.0_dp, -1.0_dp, 0.0_dp, & - 0.0_dp, 1.0_dp, 0.0_dp, & - 0.0_dp, 0.0_dp, -1.0_dp, & - 0.0_dp, 0.0_dp, 1.0_dp /),& - shape(face_normals)) - !Face normals - select case(trim(adjustl(element_type))) - case('fcc') - do i = 1, 6 - !Box face normal - face_normals(:,i) = matmul(fcc_inv, matmul(orient_inv, face_normals(:,i))) - end do - end select - !Face centroids - do i =1, 6 - - !Initialize variables - A(:) = 0.0_dp - Cx(:) = 0.0_dp - Cy = 0.0_dp - Cz = 0.0_dp - - !Calculate all terms using a projection onto the xy and xz planes and then using the 2d equation - !for centroid of a plane. This is why we calculate the x centroid twice. - do j = 1, 4 - ! A(1) = A(1) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) - - ! !Centroid in x from xy plane - ! Cx(1) = Cx(1) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) - - ! !Centroid in y from xy plane - ! Cy = Cy + (box_in_lat(2,cubic_faces(j,i))+box_in_lat(2,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) - - ! A(2) = A(2) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) - - ! !Centroid in x from xz plane - ! Cx(2) = Cx(2) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) - - ! !Centroid in z from xz plane - ! Cz = Cz + (box_in_lat(3,cubic_faces(j,i))+box_in_lat(3,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) - - ! print *, j, i, Cx, Cy, Cz, A - Cx(1) = Cx(1) + box_in_lat(1,cubic_faces(j,i))*0.25 - Cy = Cy + box_in_lat(2,cubic_faces(j,i))*0.25 - Cz = Cz + box_in_lat(3,cubic_faces(j,i))*0.25 - - end do - - ! Cx = Cx * 1/(6*A) - ! if(Cx(1) /= Cx(2)) then - ! call error_message(7) - ! end if - ! Cy = Cy* 1/(6*A(1)) - ! Cz = Cz*1/(6*A(2)) - - box_face_center(:,i) = (/ Cx(1), Cy, Cz /) - end do !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case !and one for the regular case - print *, box_bd if (esize==2) then !atomistics do iz = bd_in_lat(5)-5, bd_in_lat(6)+5 @@ -398,10 +322,141 @@ module mode_create end do end do else - continue + !If we are working with elements we have to use more complex code + + !Initialize finite element + ele(:,:) = (esize-1) * cubic_cell(:,:) + + !Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position. + !A value of true means that the coordinate is a lattice point which is within the box_bd + allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10,bd_in_lat(4)-bd_in_lat(3)+10,bd_in_lat(6)-bd_in_lat(5)+10)) + lat_points(:,:,:) = .false. + do iz = bd_in_lat(5)-5, bd_in_lat(6)+5 + do iy = bd_in_lat(3)-5, bd_in_lat(4)+5 + do ix = bd_in_lat(1)-5, bd_in_lat(2)+5 + v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) + + !Transform point back to real space for easier checking + v = matmul(orient, matmul(transform_matrix,v)) + !If within the boundaries + if(in_block_bd(v, box_bd)) then + lat_points(ix-bd_in_lat(1)+5,iy-bd_in_lat(3)+5,iz-bd_in_lat(5) + 5) = .true. + end if + end do + end do + end do + + !Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array + bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 + bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 + bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 + !Figure out where the starting point is. This is the first piont which fully contains the finite element + outerloop: do iz = 1, bd_in_array(3) + do iy = 1, bd_in_array(2) + do ix = 1, bd_in_array(1) + node_in_bd(8) = .false. + do inod = 1, 8 + vlat = ele(:,inod) + (/ ix, iy, iz /) + + !Check to see if the node resides at a position containing a lattice point within the box + if(any(vlat > shape(lat_points))) then + continue + else if(lat_points(vlat(1),vlat(2),vlat(3))) then + node_in_bd(inod) = .true. + end if + end do + + if(all(node_in_bd)) then + rzero = (/ ix, iy, iz /) + exit outerloop + end if + end do + end do + end do outerloop + + !Now build the finite element region + lat_ele_num = 0 + lat_atom_num = 0 + allocate(r_lat(3,8,numlatpoints/esize)) + + !Redefined the second 3 indices as the number of elements that fit in the bounds + do i = 1, 3 + bd_in_array(3+i) = bd_in_array(i)/esize + end do + + !Now start the element at rzero + do inod=1, 8 + ele(:,inod) = ele(:,inod) + rzero + end do + do iz = -bd_in_array(6), bd_in_array(6) + do iy = -bd_in_array(5), bd_in_array(5) + do ix = -bd_in_array(4), bd_in_array(4) + node_in_bd(:) = .false. + temp_nodes(:,:,:) = 0.0_dp + temp_lat(:,:) = 0 + do inod = 1, 8 + vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /) + !Transform point back to real space for easier checking + ! v = matmul(orient, matmul(transform_matrix,v)) + do i = 1,3 + v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5) + end do + temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) + temp_lat(:,inod) = vlat + + !Check to see if the lattice point values are greater then the array limits + if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then + continue + !If within array boundaries check to see if it is a lattice point + else if(lat_points(vlat(1),vlat(2),vlat(3))) then + node_in_bd(inod) = .true. + end if + end do + + if(all(node_in_bd)) then + lat_ele_num = lat_ele_num+1 + r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) + + !Now set all the lattice points contained within an element to false + do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:)) + do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:)) + do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:)) + lat_points(m,n,o) = .false. + end do + end do + end do + + end if + end do + end do + end do + + !Now figure out how many lattice points could not be contained in elements + print *, count(lat_points) + allocate(r_atom_lat(3,count(lat_points))) + lat_atom_num = 0 + do ix = 1, bd_in_array(3) + do iy = 1, bd_in_array(2) + do iz = 1, bd_in_array(1) + !If this point is a lattice point then save the lattice point as an atom + if (lat_points(ix,iy,iz)) then + v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) + do i = 1,3 + v(i) = v(i) + real(bd_in_lat(2*i-1) - 5) + end do + !Transform point back to real space + v = matmul(orient, matmul(transform_matrix,v)) + lat_atom_num = lat_atom_num + 1 + r_atom_lat(:,lat_atom_num) = v + end if + end do + end do + end do + + print *, lat_atom_num end if - end subroutine lattice_in_box + end subroutine build_with_rhomb subroutine error_message(errorid) From fa1cb6ce5882cbc81c8b29042ca0ff827fbd7dae Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Thu, 5 Dec 2019 08:36:22 -0500 Subject: [PATCH 4/8] Fixes to mode_create, moved basis_pos from elements to mode_create, added the mb file output style --- src/elements.f90 | 12 +++------ src/io.f90 | 62 ++++++++++++++++++++++++++++++++++++++------- src/mode_create.f90 | 22 +++++++++------- 3 files changed, 70 insertions(+), 26 deletions(-) diff --git a/src/elements.f90 b/src/elements.f90 index 21de971..3b66afb 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -30,19 +30,15 @@ module elements !Below are lattice type arrays which provide information on the general form of the elements. !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased. - + integer :: lattice_types = 0 integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type integer :: max_esize=0 !Maximum number of atoms per side of element !These variables contain information on the basis, for simplicities sake we limit !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point. !This can be easily increased with no change to efficiency - integer :: max_basisnum, basisnum(10), basis_type(10,10)!Max basis atom number, number of basis atoms in each lattice type - real(kind=dp) :: basis_pos(3,10,10) !Basis atom positions - - !Simulation cell parameters - real(kind=dp) :: box_bd(6) - + integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type + integer :: basis_type(10,10) public contains @@ -89,7 +85,6 @@ module elements max_basisnum = 0 basisnum(:) = 0 - basis_pos(:,:,:) = 0.0_dp ng_node(:) = 0 end subroutine lattice_init @@ -304,4 +299,5 @@ module elements return end subroutine rhombshape + end module elements \ No newline at end of file diff --git a/src/io.f90 b/src/io.f90 index 7850ee2..4801de1 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -3,6 +3,7 @@ module io use elements use parameters use atoms + use box implicit none @@ -56,18 +57,10 @@ module io cycle end if select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) - case('xyz') + case('xyz', 'lmp', 'vtk', 'mb') outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit - case('lmp') - outfilenum=outfilenum+1 - outfiles(outfilenum) = temp_outfile - exit - case('vtk') - outfilenum=outfilenum+1 - outfiles(outfilenum)=temp_outfile - exit case default print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", & "please input a filename with extension from following list: xyz, lmp, vtk." @@ -96,6 +89,8 @@ module io call write_lmp(outfiles(i)) case('vtk') call write_vtk(outfiles(i)) + case('mb') + call write_mb(outfiles(i)) case default print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & " is not accepted for writing. Please select from: xyz and try again" @@ -276,4 +271,53 @@ module io end do close(11) end subroutine + + subroutine write_mb(file) + + !This subroutine writes the cacmb formatted file which provides necessary information for building models + character(len=100), intent(in) :: file + + integer :: i, j, inod, ibasis + + !Open the .mb file for writing + open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + + !First write the box boundary information + !Write the global box boundaries + write(11,*) box_bd(:) + !Write the number of sub_boxes in the system + write(11,*) sub_box_num + !For every subbox write the orientation and sub box boundary + do i = 1, sub_box_num + write(11,*) sub_box_ori(:,:,i) + write(11,*) sub_box_bd(:,i) + end do + + !Write the number of atom types in the current model and all of their names + write(11,*) atom_types, (type_to_name(i), i=1, atom_types) + !Write the number of lattice_types, basisnum and number of nodes for each lattice type + write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types) + !Now for every lattice type write the basis atom types + write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types) + + !Now write the numbers of elements and atoms + write(11,*) atom_num, ele_num + + !Write out atoms first + do i = 1, atom_num + write(11,*) type_atom(i), r_atom(:,i) + end do + + !Write out the elements, this is written in two stages, one line for the element and then 1 line for + !every basis at every node + do i = 1, ele_num + write(11, *) i, lat_ele(i), size_ele(i), type_ele(i) + do inod = 1, ng_node(lat_ele(i)) + do ibasis =1, basisnum(lat_ele(i)) + write(11,*) inod, ibasis, r_node(:, ibasis, inod, i) + end do + end do + end do + + end subroutine write_mb end module io \ No newline at end of file diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 52342e2..2d41ae8 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -6,13 +6,15 @@ module mode_create use io use subroutines use elements + use box implicit none character(len=100) :: name, element_type real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), & orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3) - integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6) + integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), & + basis_pos(3,10) logical :: dup_flag, dim_flag real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) @@ -98,7 +100,7 @@ module mode_create !Add the basis atoms to the unit cell do inod = 1, max_ng_node do ibasis = 1, basisnum(1) - r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) + r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis) + origin(:) end do end do do i = 1,3 @@ -115,7 +117,7 @@ module mode_create if(lat_atom_num > 0) then do i = 1, lat_atom_num do ibasis = 1, basisnum(1) - call add_atom(basis_type(ibasis,1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis,1)) + call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) end do end do deallocate(r_atom_lat) @@ -125,7 +127,7 @@ module mode_create do i = 1, lat_ele_num do inod= 1, ng_node(1) do ibasis = 1, basisnum(1) - r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis,1) + r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis) end do end do call add_element(element_type, esize, 1, r_node_temp) @@ -258,13 +260,17 @@ module mode_create !Now normalize the orientation matrix orient = matrix_normal(orient,3) + !Set lattice_num to 1 + lattice_types = 1 + !If we haven't defined a basis then define the basis and add the default basis atom type and position if (basisnum(1) == 0) then max_basisnum = 1 basisnum(1) = 1 call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name - basis_pos(:,1,1) = 0.0_dp + basis_pos(:,1) = 0.0_dp end if + ! end subroutine subroutine build_with_rhomb(box_in_lat, transform_matrix) @@ -432,12 +438,11 @@ module mode_create end do !Now figure out how many lattice points could not be contained in elements - print *, count(lat_points) allocate(r_atom_lat(3,count(lat_points))) lat_atom_num = 0 - do ix = 1, bd_in_array(3) + do iz = 1, bd_in_array(3) do iy = 1, bd_in_array(2) - do iz = 1, bd_in_array(1) + do ix = 1, bd_in_array(1) !If this point is a lattice point then save the lattice point as an atom if (lat_points(ix,iy,iz)) then v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) @@ -453,7 +458,6 @@ module mode_create end do end do - print *, lat_atom_num end if end subroutine build_with_rhomb From 03f69c6df7942360f466b1dc7c068043ac5f314b Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Thu, 5 Dec 2019 11:03:18 -0500 Subject: [PATCH 5/8] Removed extra variables from mode_create.f90, added a new module to contain simulation box information and changed code accordingly, new grow subroutine in elements. --- src/box.f90 | 48 +++++++++++++++++++++++++++++++++++ src/elements.f90 | 61 ++++++++++++++++++++++++++++++++++++++++----- src/main.f90 | 4 ++- src/mode_create.f90 | 10 +++++--- 4 files changed, 113 insertions(+), 10 deletions(-) create mode 100644 src/box.f90 diff --git a/src/box.f90 b/src/box.f90 new file mode 100644 index 0000000..fd1586a --- /dev/null +++ b/src/box.f90 @@ -0,0 +1,48 @@ +module box + !This module contains information on the properties of the current box. + use parameters + + implicit none + + real(kind=dp) :: box_bd(6) !Global box boundaries + + !The subbox variables contain values for each subbox, being the boxes read in through some + !command. Currently only mode_merge will require sub_boxes, for mode_create it will always + !allocate to only 1 sub_box + integer :: sub_box_num = 0 + real(kind=dp), allocatable :: sub_box_ori(:,:,:) + real(kind=dp), allocatable :: sub_box_bd(:,:) + + public + contains + + subroutine box_init + !Initialize some box functions + box_bd(:) = 0.0_dp + end subroutine box_init + + subroutine alloc_sub_box(n) + !Allocate the sub_box variables + + integer, intent(in) :: n + + sub_box_num = n + allocate(sub_box_ori(3,3,n), sub_box_bd(6,n)) + end subroutine alloc_sub_box + + subroutine grow_box(temp_box_bd) + !This function takes in a temporary box boundary and adjusts the overall box boundaries + !to include it + + real(kind=dp), dimension(6), intent(in) :: temp_box_bd + + integer :: i + + do i = 1, 3 + if(temp_box_bd(2*i-1) < box_bd(2*i-1)) box_bd(2*i-1) = temp_box_bd(2*i-1) + if(temp_box_bd(2*i) > box_bd(2*i)) box_bd(2*i) = temp_box_bd(2*i) + end do + return + end subroutine grow_box + +end module box \ No newline at end of file diff --git a/src/elements.f90 b/src/elements.f90 index 3b66afb..633a50a 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -7,7 +7,6 @@ module elements implicit none !Data structures used to represent the CAC elements. Each index represents an element - integer,allocatable :: tag_ele(:) !Element tag (used to keep track of id's character(len=100), allocatable :: type_ele(:) !Element type integer, allocatable :: size_ele(:), lat_ele(:) !Element siz real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array @@ -16,7 +15,7 @@ module elements integer :: node_num=0 !Total number of nodes !Data structure used to represent atoms - integer, allocatable :: tag_atom(:), type_atom(:)!atom id + integer, allocatable :: type_atom(:)!atom type real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms @@ -119,7 +118,7 @@ module elements !Allocate element arrays if(n > 0) then - allocate(tag_ele(n), type_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), & + allocate(type_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), & stat=allostat) if(allostat > 0) then print *, "Error allocating element arrays in elements.f90 because of: ", allostat @@ -129,7 +128,7 @@ module elements if(m > 0) then !Allocate atom arrays - allocate(tag_atom(m), type_atom(m), r_atom(3,m), stat=allostat) + allocate(type_atom(m), r_atom(3,m), stat=allostat) if(allostat > 0) then print *, "Error allocating atom arrays in elements.f90 because of: ", allostat stop @@ -137,6 +136,58 @@ module elements end if end subroutine + subroutine grow_ele_arrays(n, m) + + integer, intent(in) :: n, m + + integer :: ele_size, atom_size, buffer_size + integer, allocatable :: temp_int(:) + real(kind=dp), allocatable :: temp_ele_real(:,:,:,:), temp_real(:,:) + character(len=100), allocatable :: char_temp(:) + + !The default size we grow the + buffer_size = 1024 + !Figure out the size of the atom and element arrays + ele_size = size(size_ele) + atom_size = size(type_atom) + + !Check if we need to grow the ele_size, if so grow all the variables + if ( n > size(size_ele)) then + + allocate(temp_int(n+buffer_size)) + temp_int(1:ele_size) = lat_ele + temp_int(ele_size+1:) = 0 + call move_alloc(temp_int(1:ele_size), lat_ele) + + allocate(temp_int(n+buffer_size)) + temp_int(1:ele_size) = size_ele + temp_int(ele_size+1:) = 0 + call move_alloc(temp_int(1:ele_size), size_ele) + + allocate(char_temp(n+buffer_size)) + char_temp(1:ele_size) = type_ele + call move_alloc(char_temp, type_ele) + + allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+buffer_size)) + temp_ele_real(:,:,:,1:ele_size) = r_node + temp_ele_real(:,:,:,ele_size+1:) = 0.0_dp + call move_alloc(temp_ele_real, r_node) + end if + + !Now grow atom arrays if needed + if (m > atom_size) then + allocate(temp_int(m+buffer_size)) + temp_int(1:atom_size) = type_atom + temp_int(atom_size+1:) = 0 + call move_alloc(temp_int, type_atom) + + allocate(temp_real(3,m+buffer_size)) + temp_real(:,1:atom_size) = r_atom + temp_real(:, atom_size+1:) = 0.0_dp + call move_alloc(temp_real, r_atom) + end if + end subroutine + subroutine add_element(type, size, lat, r) !Subroutine which adds an element to the element arrays integer, intent(in) :: size, lat @@ -144,7 +195,6 @@ module elements real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node) ele_num = ele_num + 1 - tag_ele(ele_num) = ele_num type_ele(ele_num) = type size_ele(ele_num) = size lat_ele(ele_num) = lat @@ -160,7 +210,6 @@ module elements real(kind=dp), intent(in), dimension(3) :: r atom_num = atom_num+1 - tag_atom(atom_num) = atom_num type_atom(atom_num) = type r_atom(:,atom_num) = r(:) diff --git a/src/main.f90 b/src/main.f90 index 7f76b6b..c7ac008 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -20,8 +20,10 @@ program main integer :: arg_num character(len=100) :: mode + !Call initialization functions call lattice_init - + call box_init + ! Command line parsing arg_num = command_argument_count() diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 2d41ae8..50b5f4c 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -135,6 +135,10 @@ module mode_create end if end if + !The last thing we do is setup the sub_box_boundaries + call alloc_sub_box(1) + sub_box_ori(:,:,1) = orient + sub_box_bd(:,1) = box_bd end subroutine create !This subroutine parses the command and pulls out information needed for mode_create subroutine parse_command() @@ -280,9 +284,9 @@ module mode_create integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space !Internal variables - integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, & - type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o - real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3) + integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), & + vlat(3), temp_lat(3,8), m, n, o + real(kind=dp) :: v(3), temp_nodes(3,1,8) real(kind=dp), allocatable :: resize_lat_array(:,:) logical, allocatable :: lat_points(:,:,:) logical :: node_in_bd(8) From d624e6ed5d1ccd1065d71f910fe8679d0f52516d Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Thu, 5 Dec 2019 11:04:49 -0500 Subject: [PATCH 6/8] Added mode_convert and documentation. Created new .mb filetype which contains all necessary information for operating on simulation cells. Created framework for reading in data files, right now only .mb style files --- README.md | 10 +++- src/Makefile | 8 +-- src/call_mode.f90 | 6 +- src/io.f90 | 133 ++++++++++++++++++++++++++++++++++++++++++- src/mode_convert.f90 | 27 +++++++++ 5 files changed, 174 insertions(+), 10 deletions(-) create mode 100644 src/mode_convert.f90 diff --git a/README.md b/README.md index 9fa35ee..34f1f31 100644 --- a/README.md +++ b/README.md @@ -91,4 +91,12 @@ origin x y z Default origin is `0 0 0`. This command just sets the origin for where the simulation cell starts building. -*Example:* `origin 10 0 1` \ No newline at end of file +*Example:* `origin 10 0 1` + +### Mode Convert + +``` +cacmb --convert infile outfile +``` + +This mode converts a file of type infile to a file of type outfile \ No newline at end of file diff --git a/src/Makefile b/src/Makefile index 58507b4..22c84eb 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,8 +1,8 @@ FC=ifort -FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone +FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin #FFLAGS=-c -mcmodel=large -Ofast -MODES=mode_create.o -OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o $(MODES) +MODES=mode_create.o mode_merge.o mode_convert.o +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o @@ -28,6 +28,6 @@ $(OBJECTS) : parameters.o atoms.o subroutines.o testfuncs.o : functions.o main.o io.o build_subroutines.o: elements.o call_mode.o : $(MODES) -$(MODES) io.o: atoms.o +$(MODES) io.o: atoms.o box.o $(MODES) main.o : io.o testfuncs.o elements.o mode_create.o: subroutines.o diff --git a/src/call_mode.f90 b/src/call_mode.f90 index 1ad4f0f..8c58d3f 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -3,6 +3,7 @@ subroutine call_mode(arg_num,mode) !mode module. use mode_create + use mode_convert use parameters implicit none @@ -12,8 +13,9 @@ subroutine call_mode(arg_num,mode) select case(mode) case('--create') - call create() - + call create + case('--convert') + call convert case default print *, "Mode ", mode, " currently not accepted. Please check documentation for ", & "accepted modes and rerun." diff --git a/src/io.f90 b/src/io.f90 index 4801de1..e5c111d 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -7,12 +7,13 @@ module io implicit none - integer :: outfilenum = 0 - character(len=100) :: outfiles(10) + integer :: outfilenum = 0, infilenum = 0 + character(len=100) :: outfiles(10), infiles(10) public contains + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutines for writing out data files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine get_out_file(filename) implicit none @@ -305,7 +306,7 @@ module io !Write out atoms first do i = 1, atom_num - write(11,*) type_atom(i), r_atom(:,i) + write(11,*) i, type_atom(i), r_atom(:,i) end do !Write out the elements, this is written in two stages, one line for the element and then 1 line for @@ -320,4 +321,130 @@ module io end do end subroutine write_mb + + !!!!!!!!!!!!! Below are subroutines for reading files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine get_in_file(filename) + + implicit none + + character(len=100), intent(in) :: filename + character(len=100) :: temp_infile + logical :: file_exists + + !If no filename is provided then this function is called with none and prompts user input + if (filename=='none') then + print *, "Please specify a filename or extension to output to:" + read(*,*) temp_infile + else + temp_infile = filename + end if + + !Infinite loop which only exists if user provides valid filetype + do while(.true.) + + !Check to see if file exists, if it does then ask user if they would like to overwrite the file + inquire(file=trim(temp_infile), exist=file_exists) + if (.not.file_exists) then + print *, "The file ", filename, " does not exist. Please input a filename that exists" + read(*,*) temp_infile + cycle + end if + + select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) + case('xyz', 'lmp', 'vtk', 'mb') + infilenum=infilenum+1 + infiles(infilenum) = temp_infile + exit + case default + print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), "not currently accepted. ", & + "please input a filename with extension from following list: mb." + read(*,*) temp_infile + + end select + end do + + end subroutine get_in_file + + subroutine read_in + !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine + + integer :: i + + do i = 1, infilenum + !Pull out the extension of the file and call the correct write subroutine + select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:)))) + case('mb') + call read_mb(infiles(i)) + case default + print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & + " is not accepted for writing. Please select from: mb and try again" + stop + + end select + end do + end subroutine read_in + + subroutine read_mb(file) + !This subroutine reads in an mb file for operation + + character(len=100), intent(in) :: file + + integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles + character(len=100) :: etype + real(kind=dp) :: temp_box_bd(6), r(3) + real(kind=dp), allocatable :: r_innode(:,:,:) + !First open the file + open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') + + !Read in the box boundary and grow the current active box bd + read(11, *) temp_box_bd(:) + call grow_box(temp_box_bd) + + !Read in the number of sub_boxes and allocate the variables + read(11, *) n + call alloc_sub_box(n) + + !Read in subbox orientations and boundaries + do i = 1, sub_box_num + !Read in orientation with column major ordering + read(11,*) ((sub_box_ori(j, k, i), j = 1, 3), k = 1, 3) + !Read in subbox boundaries + read(11,*) sub_box_bd(:,i) + end do + + !Read in the number of atom types and all their names + read(11, *) atom_types, (type_to_name(i), i = 1, atom_types) + !Read the number of lattice types, basisnum, and number of nodes for each lattice type + read(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types) + !Define max_ng_node and max_basis_num + max_basisnum = maxval(basisnum) + max_ng_node = maxval(ng_node) + !Read the basis atom types for every lattice + read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types) + + !Read number of elements and atoms and allocate arrays + read(11, *) in_atoms,in_eles + call grow_ele_arrays(in_eles, in_atoms) + allocate(r_innode(3,max_basisnum, max_ng_node)) + + !Read the atoms + do i = 1, in_atoms + read(11,*) j, type, r(:) + call add_atom(type, r) + end do + + !Read the elements + do i = 1, in_eles + read(11, *) n, type, size, etype + do inod = 1, ng_node(type) + do ibasis =1, basisnum(type) + read(11,*) j, k, r_innode(:, ibasis, inod) + end do + end do + + call add_element(etype, size, type, r_innode) + end do + + end subroutine read_mb end module io \ No newline at end of file diff --git a/src/mode_convert.f90 b/src/mode_convert.f90 new file mode 100644 index 0000000..b46b07f --- /dev/null +++ b/src/mode_convert.f90 @@ -0,0 +1,27 @@ +module mode_convert + + use parameters + use box + use elements + use io + + public + contains + + subroutine convert + !This subroutine converts a single input file from one format to another + character(len=100) :: infile, outfile + + !We have to allocate the element and atom arrays with a size of 1 for the read in code to work + call alloc_ele_arrays(1,1) + !First read in the file + call get_command_argument(2, infile) + call get_in_file(infile) + call read_in + + !Now get the outfile, writing is done after all the codes complete + call get_command_argument(3, outfile) + call get_out_file(outfile) + + end subroutine convert +end module mode_convert \ No newline at end of file From fb2abc60d1f9d2d2467931e784190903a92ae299 Mon Sep 17 00:00:00 2001 From: Alex Date: Sat, 7 Dec 2019 13:38:52 -0500 Subject: [PATCH 7/8] Updates to box variables adding new sub_box variables and propagating changes through all of the modes --- README.md | 14 +++++++++++++- src/box.f90 | 35 +++++++++++++++++++++++++++++++---- src/elements.f90 | 16 ++++++++-------- src/mode_convert.f90 | 5 +++-- src/mode_create.f90 | 4 ++++ 5 files changed, 59 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index 34f1f31..d390177 100644 --- a/README.md +++ b/README.md @@ -99,4 +99,16 @@ Default origin is `0 0 0`. This command just sets the origin for where the simul cacmb --convert infile outfile ``` -This mode converts a file of type infile to a file of type outfile \ No newline at end of file +This mode converts a file `infile` to a file of `outfile`. The extensions determine the conversion process. + +### Mode Merge + +``` +cacmb --merge dim N infiles outfile +``` + +This mode merges multiple data files and creates one big simulation cell. The parameters are: + +`N` - The number of files which are being read + +`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command. \ No newline at end of file diff --git a/src/box.f90 b/src/box.f90 index fd1586a..ad5fe55 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -10,8 +10,9 @@ module box !command. Currently only mode_merge will require sub_boxes, for mode_create it will always !allocate to only 1 sub_box integer :: sub_box_num = 0 - real(kind=dp), allocatable :: sub_box_ori(:,:,:) - real(kind=dp), allocatable :: sub_box_bd(:,:) + integer, allocatable :: sub_box_array_bd(:,:,:)!Boundaries in the atom and element arrays for each sub_box + real(kind=dp), allocatable :: sub_box_ori(:,:,:)!Orientations for each of the subboxes + real(kind=dp), allocatable :: sub_box_bd(:,:)!Boundaries for each of the sub_boxes public contains @@ -26,10 +27,36 @@ module box integer, intent(in) :: n - sub_box_num = n - allocate(sub_box_ori(3,3,n), sub_box_bd(6,n)) + allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n)) + end subroutine alloc_sub_box + subroutine grow_sub_box(n) + !Grows sub box arrays, this is only called when a new file is read in + integer, intent(in) :: n + + integer, allocatable :: temp_array_bd(:,:,:), temp_file(:) + real(kind=dp), allocatable :: temp_ori(:,:,:), temp_bd(:,:) + !Allocate temporary arrays + allocate(temp_ori(3,3,sub_box_num+n),temp_bd(6,sub_box_num+n), & + temp_array_bd(2,2,sub_box_num+n), temp_file(sub_box_num+n)) + + !Move allocation for all sub_box_arrays + temp_ori(:,:,1:sub_box_num) = sub_box_ori + temp_ori(:,:,sub_box_num+1:) = 0.0_dp + call move_alloc(temp_ori, sub_box_ori) + + temp_bd(:, 1:sub_box_num) = sub_box_bd + temp_bd(:, sub_box_num+1:) = 0.0_dp + call move_alloc(temp_bd, sub_box_bd) + + temp_array_bd(:,:,1:sub_box_num) = sub_box_array_bd + temp_array_bd(:,:,sub_box_num+1:) = 0.0_dp + call move_alloc(temp_array_bd, sub_box_array_bd) + + return + end subroutine grow_sub_box + subroutine grow_box(temp_box_bd) !This function takes in a temporary box boundary and adjusts the overall box boundaries !to include it diff --git a/src/elements.f90 b/src/elements.f90 index 633a50a..28e0e5e 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -152,36 +152,36 @@ module elements atom_size = size(type_atom) !Check if we need to grow the ele_size, if so grow all the variables - if ( n > size(size_ele)) then + if ( n+ele_num > size(size_ele)) then - allocate(temp_int(n+buffer_size)) + allocate(temp_int(n+ele_num+buffer_size)) temp_int(1:ele_size) = lat_ele temp_int(ele_size+1:) = 0 call move_alloc(temp_int(1:ele_size), lat_ele) - allocate(temp_int(n+buffer_size)) + allocate(temp_int(n+ele_num+buffer_size)) temp_int(1:ele_size) = size_ele temp_int(ele_size+1:) = 0 call move_alloc(temp_int(1:ele_size), size_ele) - allocate(char_temp(n+buffer_size)) + allocate(char_temp(n+ele_num+buffer_size)) char_temp(1:ele_size) = type_ele call move_alloc(char_temp, type_ele) - allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+buffer_size)) + allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_num+buffer_size)) temp_ele_real(:,:,:,1:ele_size) = r_node temp_ele_real(:,:,:,ele_size+1:) = 0.0_dp call move_alloc(temp_ele_real, r_node) end if !Now grow atom arrays if needed - if (m > atom_size) then - allocate(temp_int(m+buffer_size)) + if (m+atom_num > atom_size) then + allocate(temp_int(m+atom_num+buffer_size)) temp_int(1:atom_size) = type_atom temp_int(atom_size+1:) = 0 call move_alloc(temp_int, type_atom) - allocate(temp_real(3,m+buffer_size)) + allocate(temp_real(3,m+atom_num+buffer_size)) temp_real(:,1:atom_size) = r_atom temp_real(:, atom_size+1:) = 0.0_dp call move_alloc(temp_real, r_atom) diff --git a/src/mode_convert.f90 b/src/mode_convert.f90 index b46b07f..d6e1d9b 100644 --- a/src/mode_convert.f90 +++ b/src/mode_convert.f90 @@ -11,13 +11,14 @@ module mode_convert subroutine convert !This subroutine converts a single input file from one format to another character(len=100) :: infile, outfile - + real(kind = dp) :: temp_box_bd(6) !We have to allocate the element and atom arrays with a size of 1 for the read in code to work call alloc_ele_arrays(1,1) !First read in the file call get_command_argument(2, infile) call get_in_file(infile) - call read_in + call read_in(1, (/0.0_dp,0.0_dp,0.0_dp/), temp_box_bd) + call grow_box(temp_box_bd) !Now get the outfile, writing is done after all the codes complete call get_command_argument(3, outfile) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 50b5f4c..b82cc82 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -137,8 +137,12 @@ module mode_create !The last thing we do is setup the sub_box_boundaries call alloc_sub_box(1) + sub_box_num = 1 sub_box_ori(:,:,1) = orient sub_box_bd(:,1) = box_bd + sub_box_array_bd(1,:,1) = 1 + sub_box_array_bd(2,1,1) = atom_num + sub_box_array_bd(2,2,1) = ele_num end subroutine create !This subroutine parses the command and pulls out information needed for mode_create subroutine parse_command() From 3c7461f094c120de40b48a84ec3115d7e8481ca2 Mon Sep 17 00:00:00 2001 From: Alex Date: Sat, 7 Dec 2019 13:39:29 -0500 Subject: [PATCH 8/8] Added mode merge, adjusted how file reading works to accomodate model merge --- src/call_mode.f90 | 7 +++- src/io.f90 | 72 +++++++++++++++++++++----------- src/mode_merge.f90 | 102 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 154 insertions(+), 27 deletions(-) create mode 100644 src/mode_merge.f90 diff --git a/src/call_mode.f90 b/src/call_mode.f90 index 8c58d3f..f571520 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -4,6 +4,7 @@ subroutine call_mode(arg_num,mode) use mode_create use mode_convert + use mode_merge use parameters implicit none @@ -15,9 +16,11 @@ subroutine call_mode(arg_num,mode) case('--create') call create case('--convert') - call convert + call convert + case('--merge') + call merge case default - print *, "Mode ", mode, " currently not accepted. Please check documentation for ", & + print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", & "accepted modes and rerun." stop 3 diff --git a/src/io.f90 b/src/io.f90 index e5c111d..20fac52 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -278,7 +278,7 @@ module io !This subroutine writes the cacmb formatted file which provides necessary information for building models character(len=100), intent(in) :: file - integer :: i, j, inod, ibasis + integer :: i, j, k, inod, ibasis !Open the .mb file for writing open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') @@ -288,10 +288,11 @@ module io write(11,*) box_bd(:) !Write the number of sub_boxes in the system write(11,*) sub_box_num - !For every subbox write the orientation and sub box boundary + !For every subbox write the orientation, sub box boundary, and sub_box_array_bds do i = 1, sub_box_num write(11,*) sub_box_ori(:,:,i) write(11,*) sub_box_bd(:,i) + write(11,*) ((sub_box_array_bd(j,k,i), j = 1, 2), k = 1, 2) end do !Write the number of atom types in the current model and all of their names @@ -346,7 +347,7 @@ module io !Check to see if file exists, if it does then ask user if they would like to overwrite the file inquire(file=trim(temp_infile), exist=file_exists) if (.not.file_exists) then - print *, "The file ", filename, " does not exist. Please input a filename that exists" + print *, "The file ", trim(adjustl(filename)), " does not exist. Please input a filename that exists" read(*,*) temp_infile cycle end if @@ -366,52 +367,70 @@ module io end subroutine get_in_file - subroutine read_in + subroutine read_in(i, displace, temp_box_bd) !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine - integer :: i + integer, intent(in) :: i + real(kind=dp), dimension(3), intent(in) :: displace + real(kind=dp), dimension(6), intent(out) :: temp_box_bd - do i = 1, infilenum - !Pull out the extension of the file and call the correct write subroutine - select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:)))) - case('mb') - call read_mb(infiles(i)) - case default - print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & - " is not accepted for writing. Please select from: mb and try again" - stop + !Pull out the extension of the file and call the correct write subroutine + select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:)))) + case('mb') + call read_mb(infiles(i), displace, temp_box_bd) + case default + print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & + " is not accepted for writing. Please select from: mb and try again" + stop + + end select - end select - end do end subroutine read_in - subroutine read_mb(file) + subroutine read_mb(file, displace, temp_box_bd) !This subroutine reads in an mb file for operation character(len=100), intent(in) :: file + real(kind=dp), dimension(3), intent(in) :: displace + real(kind = dp), dimension(6), intent(out) :: temp_box_bd integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles character(len=100) :: etype - real(kind=dp) :: temp_box_bd(6), r(3) + real(kind=dp) :: r(3), newdisplace(3) real(kind=dp), allocatable :: r_innode(:,:,:) !First open the file open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') !Read in the box boundary and grow the current active box bd read(11, *) temp_box_bd(:) - call grow_box(temp_box_bd) + do i = 1, 3 + newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i) + temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i) + end do + !Read in the number of sub_boxes and allocate the variables read(11, *) n - call alloc_sub_box(n) + + if (sub_box_num == 0) then + call alloc_sub_box(n) + else + call grow_sub_box(n) + end if !Read in subbox orientations and boundaries - do i = 1, sub_box_num + do i = 1, n !Read in orientation with column major ordering - read(11,*) ((sub_box_ori(j, k, i), j = 1, 3), k = 1, 3) + read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 3), k = 1, 3) !Read in subbox boundaries - read(11,*) sub_box_bd(:,i) + read(11,*) sub_box_bd(:,sub_box_num+i) + sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:) + !Read in sub_box_array_bd + read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 2), k = 1, 2) + end do + sub_box_num = sub_box_num + n !Read in the number of atom types and all their names read(11, *) atom_types, (type_to_name(i), i = 1, atom_types) @@ -424,14 +443,14 @@ module io read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types) !Read number of elements and atoms and allocate arrays - read(11, *) in_atoms,in_eles + read(11, *) in_atoms, in_eles call grow_ele_arrays(in_eles, in_atoms) allocate(r_innode(3,max_basisnum, max_ng_node)) !Read the atoms do i = 1, in_atoms read(11,*) j, type, r(:) - call add_atom(type, r) + call add_atom(type, r+newdisplace) end do !Read the elements @@ -440,11 +459,14 @@ module io do inod = 1, ng_node(type) do ibasis =1, basisnum(type) read(11,*) j, k, r_innode(:, ibasis, inod) + r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace end do end do call add_element(etype, size, type, r_innode) end do + !Close the file being read + close(11) end subroutine read_mb end module io \ No newline at end of file diff --git a/src/mode_merge.f90 b/src/mode_merge.f90 new file mode 100644 index 0000000..0919b70 --- /dev/null +++ b/src/mode_merge.f90 @@ -0,0 +1,102 @@ +module mode_merge + !This module contains the code needed for merging several .mb files together + + use parameters + use atoms + use io + use subroutines + use elements + + character(len=4) :: dim + integer :: in_num + + public + contains + subroutine merge + + integer :: i + real(kind=dp) :: displace(3), temp_box_bd(6) + + !First we parse the merge command + call parse_command + + !Now loop over all files and stack them + do i = 1, in_num + displace(:) = 0.0_dp + if ((i==1).or.(trim(adjustl(dim)) == 'none')) then + call read_in(i, displace, temp_box_bd) + call grow_box(temp_box_bd) + else + select case(trim(adjustl(dim))) + case('x') + displace(1) = box_bd(2) + case('y') + displace(2) = box_bd(4) + case('z') + displace(3) = box_bd(6) + end select + + call read_in(i, displace, temp_box_bd) + call grow_box(temp_box_bd) + end if + end do + + return + end subroutine merge + + subroutine parse_command + + character(len=100) :: textholder + integer :: i, stat, arglen, arg_pos + + !Get dimension to concatenate along + call get_command_argument(2, dim, arglen) + if (arglen == 0) STOP "Missing dim in mode_merge command" + select case(trim(adjustl(dim))) + case('x','y','z','none') + continue + case default + print *, dim, " not accepted as a dimension in mode_merge" + stop 3 + end select + !Get number of files to read in + call get_command_argument(3, textholder, arglen) + if (arglen == 0) STOP "Number of files to read missing in mode_merge command" + read(textholder, *, iostat = stat) in_num + if (stat > 0) STOP "Error reading number of files in, must be integer" + + !Now loop and pull out all files + do i = 1, in_num + call get_command_argument(3+i, textholder, arglen) + if (arglen == 0) STOP "Fewer files to read in then specified" + !Make sure this file is readable + call get_in_file(textholder) + end do + + !Set argpos accordingly + arg_pos = 3+in_num+1 + !Now options loop + !Check for optional keywords + do while(.true.) + if(arg_pos > command_argument_count()) exit + !Pull out the next argument which should either be a keyword or an option + call get_command_argument(arg_pos, textholder) + textholder=adjustl(textholder) + arg_pos=arg_pos+1 + + !Choose what to based on what the option string is + select case(trim(textholder)) + + case default + !Check to see if it is an option command, if so then mode_create must be finished + if(textholder(1:1) == '-') then + exit + !Check to see if a filename was passed + elseif(scan(textholder,'.',.true.) > 0) then + call get_out_file(textholder) + end if + end select + end do + + end subroutine parse_command +end module mode_merge \ No newline at end of file