Latest working version of cacmb

development
Alex Selimov 3 years ago
parent fbaca5859b
commit f0fd76f12d

@ -1162,6 +1162,267 @@ END SELECT
!
END SUBROUTINE ATOMMASSSPECIES
!
subroutine realatomspecies(smass,species)
!
IMPLICIT NONE
CHARACTER(LEN=2),INTENT(OUT):: species
INTEGER:: mass
REAL(dp),INTENT(IN):: smass
if(mass_is_equal(smass,1.008d0)) then
species = 'H'
else if(mass_is_equal(smass,2.014101777d0)) then
species = 'D'
else if(mass_is_equal(smass,4.002602d0)) then
species = 'He'
else if(mass_is_equal(smass,6.94d0)) then
species = 'Li'
else if(mass_is_equal(smass,9.012182d0)) then
species = 'Be'
else if(mass_is_equal(smass,10.81d0)) then
species = 'B'
else if(mass_is_equal(smass,12.011d0)) then
species='C'
else if(mass_is_equal(smass,14.007d0)) then
species='N'
else if(mass_is_equal(smass,15.999d0)) then
species='O'
else if(mass_is_equal(smass,18.9984032d0)) then
species='F'
else if(mass_is_equal(smass,20.1797d0)) then
species='Ne'
!
! n=3
else if(mass_is_equal(smass,22.98976928d0)) then
species = 'Na'
else if(mass_is_equal(smass,24.305d0)) then
species='Mg'
else if(mass_is_equal(smass,26.9815386d0)) then
species='Al'
else if(mass_is_equal(smass,28.085d0)) then
species='Si'
else if(mass_is_equal(smass,30.973762d0)) then
species='P'
else if(mass_is_equal(smass,32.06d0)) then
species='S'
else if(mass_is_equal(smass,35.45d0)) then
species='Cl'
else if(mass_is_equal(smass,39.948d0)) then
species='Ar'
!
! n=4
else if(mass_is_equal(smass,39.0983d0)) then
species='K'
else if(mass_is_equal(smass,40.078d0)) then
species='Ca'
else if(mass_is_equal(smass,44.955912d0)) then
species='Sc'
else if(mass_is_equal(smass,47.867d0)) then
species='Ti'
else if(mass_is_equal(smass,50.9415d0)) then
species='V'
else if(mass_is_equal(smass,51.9961d0)) then
species='Cr'
else if(mass_is_equal(smass,54.938045d0)) then
species='Mn'
else if(mass_is_equal(smass,55.845d0)) then
species='Fe'
else if(mass_is_equal(smass,58.933195d0)) then
species='Co'
else if(mass_is_equal(smass,58.6934d0)) then
species='Ni'
else if(mass_is_equal(smass,63.546d0)) then
species='Cu'
else if(mass_is_equal(smass,65.38d0)) then
species='Zn'
else if(mass_is_equal(smass,69.723d0)) then
species='Ga'
else if(mass_is_equal(smass,72.63d0)) then
species='Ge'
else if(mass_is_equal(smass,74.9216d0)) then
species='As'
else if(mass_is_equal(smass,78.96d0)) then
species='Se'
else if(mass_is_equal(smass,79.904d0)) then
species='Br'
else if(mass_is_equal(smass,83.798d0)) then
species='Kr'
!
! n=5
else if(mass_is_equal(smass,85.4678d0)) then
species='Rb'
else if(mass_is_equal(smass,87.62d0)) then
species='Sr'
else if(mass_is_equal(smass,88.90585d0)) then
species='Y'
else if(mass_is_equal(smass,91.224d0)) then
species='Zr'
else if(mass_is_equal(smass,92.90638d0)) then
species='Nb'
else if(mass_is_equal(smass,95.96d0)) then
species='Mo'
else if(mass_is_equal(smass,98.906d0)) then
species='Tc'
else if(mass_is_equal(smass,101.07d0)) then
species='Ru'
else if(mass_is_equal(smass,102.90550d0)) then
species='Rh'
else if(mass_is_equal(smass,106.42d0)) then
species='Pd'
else if(mass_is_equal(smass,107.8682d0)) then
species='Ag'
else if(mass_is_equal(smass,112.411d0)) then
species='Cd'
else if(mass_is_equal(smass,114.818d0)) then
species='In'
else if(mass_is_equal(smass,118.71d0)) then
species='Sn'
else if(mass_is_equal(smass,121.76d0)) then
species='Sb'
else if(mass_is_equal(smass,127.60d0)) then
species='Te'
else if(mass_is_equal(smass,126.90447d0)) then
species='I'
else if(mass_is_equal(smass,131.293d0)) then
species='Xe'
!
! n=6
else if(mass_is_equal(smass,132.9054519d0)) then
species='Cs'
else if(mass_is_equal(smass,137.327d0)) then
species='Ba'
! Lanthanides
else if(mass_is_equal(smass,138.90547d0)) then
species='La'
else if(mass_is_equal(smass,140.116d0)) then
species='Ce'
else if(mass_is_equal(smass,140.90765d0)) then
species='Pr'
else if(mass_is_equal(smass,144.242d0)) then
species='Nd'
else if(mass_is_equal(smass,144.91d0)) then
species='Pm'
else if(mass_is_equal(smass,150.36d0)) then
species='Sm'
else if(mass_is_equal(smass,151.964d0)) then
species='Eu'
else if(mass_is_equal(smass,157.25d0)) then
species='Gd'
else if(mass_is_equal(smass,158.92535d0)) then
species='Tb'
else if(mass_is_equal(smass,162.50d0)) then
species='Dy'
else if(mass_is_equal(smass,164.93032d0)) then
species='Ho'
else if(mass_is_equal(smass,167.259d0)) then
species='Er'
else if(mass_is_equal(smass,168.93421d0)) then
species='Tm'
else if(mass_is_equal(smass,173.054d0)) then
species='Yb'
else if(mass_is_equal(smass,174.9668d0)) then
species='Lu'
! End of Lanthanides
else if(mass_is_equal(smass,178.49d0)) then
species='Hf'
else if(mass_is_equal(smass,180.94788d0)) then
species='Ta'
else if(mass_is_equal(smass,183.84d0)) then
species='W'
else if(mass_is_equal(smass,186.207d0)) then
species='Re'
else if(mass_is_equal(smass,190.23d0)) then
species='Os'
else if(mass_is_equal(smass,192.217d0)) then
species='Ir'
else if(mass_is_equal(smass,195.084d0)) then
species='Pt'
else if(mass_is_equal(smass,196.966569d0)) then
species='Au'
else if(mass_is_equal(smass,200.59d0)) then
species='Hg'
else if(mass_is_equal(smass,204.38d0)) then
species='Tl'
else if(mass_is_equal(smass,207.2d0)) then
species='Pb'
else if(mass_is_equal(smass,208.9804d0)) then
species='Bi'
else if(mass_is_equal(smass,209.98d0)) then
species='Po'
else if(mass_is_equal(smass,209.99d0)) then
species='At'
else if(mass_is_equal(smass,222.02d0)) then
species='Rn'
!
! n=7
else if(mass_is_equal(smass,233.d0)) then
species='Fr'
else if(mass_is_equal(smass,226.d0)) then
species='Ra'
! Actinides
else if(mass_is_equal(smass,227.d0)) then
species='Ac'
else if(mass_is_equal(smass,232.03806d0)) then
species='Th'
else if(mass_is_equal(smass,231.03588d0)) then
species='Pa'
else if(mass_is_equal(smass,238.02891d0)) then
species='U'
else if(mass_is_equal(smass,237.d0)) then
species='Np'
else if(mass_is_equal(smass,244.d0)) then
species='Pu'
else if(mass_is_equal(smass,243.d0)) then
species='Am'
else if(mass_is_equal(smass,247.d0)) then
species='Cm'
else if(mass_is_equal(smass,247.d0)) then
species='Bk'
else if(mass_is_equal(smass,251.d0)) then
species='Cf'
else if(mass_is_equal(smass,252.d0)) then
species='Es'
else if(mass_is_equal(smass,257.d0)) then
species='Fm'
else if(mass_is_equal(smass,258.d0)) then
species='Md'
else if(mass_is_equal(smass,259.d0)) then
species='No'
else if(mass_is_equal(smass,262.d0)) then
species='Lr'
! End of actinides
else if(mass_is_equal(smass,265.d0)) then
species='Rf'
else if(mass_is_equal(smass,268.d0)) then
species='Db'
else if(mass_is_equal(smass,271.d0)) then
species='Sg'
else if(mass_is_equal(smass,270.d0)) then
species='Bh'
else if(mass_is_equal(smass,277.d0)) then
species='Hs'
else if(mass_is_equal(smass,276.d0)) then
species='Mt'
else if(mass_is_equal(smass,281.d0)) then
species='Ds'
else if(mass_is_equal(smass,280.d0)) then
species='Rg'
else if(mass_is_equal(smass,285.17d0)) then
species='Cn'
else if(mass_is_equal(smass,284.d0)) then
species='Uu'
else if(mass_is_equal(smass,289.d0)) then
species='Fl'
else if(mass_is_equal(smass,288.d0)) then
species='Mc'
else if(mass_is_equal(smass,293.d0)) then
species='Lv'
else if(mass_is_equal(smass,294.d0)) then
species='Ts'
else if(mass_is_equal(smass,294.d0)) then
species='Og'
!
end if
end subroutine realatomspecies
END MODULE atoms

@ -14,7 +14,7 @@ module elements
integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size
real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
!Element result data structures
real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:)
real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:), vel_node(:,:,:,:)
integer, save :: ele_num !Number of elements
integer, save :: node_num !Total number of nodes
@ -26,7 +26,7 @@ module elements
real(kind =dp),allocatable :: r_atom(:,:) !atom position
integer :: atom_num=0 !Number of atoms
!Atom result data structures information
real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:)
real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:), vel_atom(:,:)
!Mapping atom type to provided name
character(len=2), dimension(10) :: type_to_name
@ -42,7 +42,7 @@ module elements
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
real(kind=dp) :: lapa(10)
real(kind=dp) :: lapa(10), masses(10)
!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.
@ -919,7 +919,6 @@ do i = 1, atom_num
integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays
integer :: allostat
print *, max_ng_node
!Allocate element arrays
if (n > 0) then
if (allocated(force_node)) then
@ -952,6 +951,38 @@ do i = 1, atom_num
end subroutine
subroutine alloc_vel_arrays(n,m)
!This subroutine used to provide initial allocation for the atom and element data 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
if (allocated(vel_node)) then
deallocate(vel_node)
end if
allocate(vel_node(3,max_basisnum,max_ng_node,n), stat=allostat)
if(allostat > 0) then
print *, "Error allocating element data arrays in mode_metric because of:", allostat
stop
end if
end if
if (m > 0) then
if (allocated(vel_atom)) then
deallocate(vel_atom)
end if
allocate(vel_atom(3,m), stat=allostat)
if(allostat > 0) then
print *, "Error allocating atom data arrays in mode_metric because of:", allostat
stop
end if
end if
end subroutine alloc_vel_arrays
subroutine add_atom_data(ia, eng, force, virial)
!Function which sets the atom data arrays
integer, intent(in) :: ia

@ -254,6 +254,19 @@ END FUNCTION StrDnCase
return
end function is_equal
function mass_is_equal(A, B)
!This function checks if too numbers are equal within a tolerance
real(kind=dp), intent(in) :: A, B
logical :: mass_is_equal
if((A>(B - 10.0_dp**(-2))).and.(A < (B+10.0_dp**(-2)))) then
mass_is_equal = .true.
else
mass_is_equal = .false.
end if
return
end function mass_is_equal
pure function unitvec(n,vec)
integer, intent(in) :: n
real(kind=dp), intent(in) :: vec(n)

@ -520,6 +520,8 @@ module io
12 format('coarse-grained domain')
13 format('ie basis_num ng_node esize'/ 'ip ibasis type x y z')
14 format('atomistic domain' / 'ia type_atom x y z')
15 format('ie basis_num ng_node esize'/ 'ip ibasis type x y z velx vely velz')
16 format('atomistic domain' / 'ia type_atom x y z velx vely velz')
19 format('max nodes per element and basis per nodes' / 2i16)
20 format('max interpo per element' / i16)
21 format('atom types to elements')
@ -562,7 +564,11 @@ module io
if(ele_num > 0) then
write(11,12)
ip = 0
if(allocated(vel_node)) then
write(11,15)
else
write(11,13)
end if
do i = 1, ele_num
select case(type_ele(i))
case('fcc','bcc')
@ -570,22 +576,39 @@ module io
case('20fcc')
write(11, '(4i16)') i, basisnum(lat_ele(i)), 3, (size_ele(i)-1)
end select
if(allocated(vel_node)) then
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
ip = ip + 1
write(11, '(3i16, 6f23.15)') ip, ibasis, basis_type(ibasis, lat_ele(i)), r_node(:, ibasis, inod, i), &
vel_node(:,ibasis,inod,i)
end do
end do
else
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
ip = ip + 1
write(11, '(3i16, 3f23.15)') ip, ibasis, basis_type(ibasis, lat_ele(i)), r_node(:, ibasis, inod, i)
end do
end do
end if
end do
end if
!Now write the atomic information
if(atom_num /= 0) then
if(allocated(vel_atom)) then
write(11, 16)
do i = 1, atom_num
write(11, '(2i16, 6f23.15)') i, type_atom(i), r_atom(:,i), vel_atom(:, i)
end do
else
write(11,14)
do i = 1, atom_num
write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i)
end do
end if
end if
close(11)
end subroutine write_pycac
@ -730,9 +753,10 @@ module io
real(kind = dp), dimension(6), intent(out) :: temp_box_bd
integer :: i, j, k, l, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
new_type_to_type(10), new_lattice_types, sbox, new_lattice_map(10)
new_type_to_type(10), new_lattice_types, sbox, new_lattice_map(10), btype_map(10,10), bnum_map(10), &
node_map(10)
character(len=100) :: etype
real(kind=dp) :: r(3), newdisplace(3)
real(kind=dp) :: r(3), newdisplace(3), lapa_map(10)
real(kind=dp), allocatable :: r_innode(:,:,:)
character(len = 2) :: new_type_to_name(10)
!First open the file
@ -782,45 +806,29 @@ module io
call add_atom_type(new_type_to_name(i), new_type_to_type(i), all_new)
end do
!Read the number of lattice types, basisnum, and number of nodes for each lattice type
read(11,*) new_lattice_types, (basisnum(i), i = lattice_types+1, lattice_types+new_lattice_types), &
(ng_node(i), i = lattice_types+1, lattice_types+new_lattice_types)
!Define max_ng_node and max_basis_num
max_basisnum = maxval(basisnum)
max_ng_node = maxval(ng_node)
read(11,*) new_lattice_types, (bnum_map(i), i = 1, new_lattice_types),(node_map(i), i =1, new_lattice_types)
!Read the basis atom types for every lattice
read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = lattice_types+1, lattice_types+new_lattice_types)
read(11,*) ((btype_map(i,j), i = 1, bnum_map(j)), j = 1, new_lattice_types)
!Convert the basis_atom types
do j = lattice_types+1, lattice_types+new_lattice_types
do i = 1, basisnum(j)
basis_type(i,j) = new_type_to_type(basis_type(i,j))
do j =1, new_lattice_types
do i = 1, bnum_map(j)
btype_map(i,j) = new_type_to_type(btype_map(i,j))
end do
end do
!Read the lattice parameters for every lattice type
read(11,*) (lapa(i), i = lattice_types+1, lattice_types+new_lattice_types)
read(11,*) (lapa_map(i), i = 1, new_lattice_types)
!Now we loop over all new lattice types and check to see if they are exactly the same as any old lattice types
!If they are then we map the new type to the old type.
k = lattice_types + 1
new_lattice_map(:) = 0
new_loop:do i = 1, new_lattice_types
old_loop:do j = 1, lattice_types
!First check all the lattice level variables
if ((basisnum(lattice_types + i) == basisnum(j)).and. &
(ng_node(lattice_types + i) == ng_node(j))) then
!Now check the basis level variables
do ibasis =1, basisnum(i)
if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then
cycle old_loop
end if
end do
new_lattice_map(i) = j
cycle new_loop
end if
end do old_loop
new_lattice_map(i) = k
k = k+1
call lattice_map(bnum_map(i), btype_map(:,i), node_map(i), lapa_map(i), new_lattice_map(i))
end do new_loop
!Define max_ng_node and max_basis_num
max_basisnum = maxval(basisnum)
max_ng_node = maxval(ng_node)
!Read number of elements and atoms and allocate arrays
read(11, *) in_atoms, in_eles
call grow_ele_arrays(in_eles, in_atoms)
@ -909,7 +917,7 @@ module io
!Read define atom_types by mass
print *, j
do i = 1, j
call atommassspecies(atomic_masses(i), atomic_element)
call realatomspecies(atomic_masses(i), atomic_element)
call add_atom_type(atomic_element, atom_type_map(i), all_new)
print *, i, atom_type_map(i)
end do
@ -1009,13 +1017,15 @@ module io
!Internal Variables
integer :: i, j, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, &
id, type_node, ilat, esize, tag, type, bnum, n, ibasis, ip, atom_type_map(100)
real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), &
ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), atomic_masses(10)
real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), vela(3),&
ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), vele(3,10,20), atomic_masses(10)
character(len=100) :: textholder, fcc
character(len=1000) :: line
character(len=2) :: atomic_element
logical :: read_vel
read_vel =.false.
open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
!Now initialize some important variables if they aren't defined
@ -1063,7 +1073,7 @@ module io
!Read define atom_types by mass
do i = 1, j
call atommassspecies(atomic_masses(i), atomic_element)
call realatomspecies(atomic_masses(i), atomic_element)
call add_atom_type(atomic_element, atom_type_map(i), all_new)
end do
@ -1077,12 +1087,22 @@ module io
if(in_atoms > 0 ) then
!Read atom header
read(11,*) textholder
read(11,'(a)') line
read(line, *) (textholder, i=1, 18)
if(textholder=="velx") then
call alloc_vel_arrays(in_eles, in_atoms)
read_vel = .true.
end if
do ia = 1, in_atoms
read(11,'(a)') line(:)
if(read_vel) then
read(line,*) tag, type, ra(:), ea, fa(:), va(:), vela(:)
else
read(line,*) tag, type, ra(:), ea, fa(:), va(:)
end if
call add_atom(tag, atom_type_map(type), sub_box_num, ra)
call add_atom_data(atom_num, ea, fa, va)
if(read_vel) vel_atom(:,atom_num) = vela
end do
end if
@ -1096,9 +1116,16 @@ module io
do ie =1, in_eles
read(11,*) tag, n, bnum, esize
inbtypes(:) = 0
if(read_vel) then
do inod =1, n*bnum
read(11,*) ip, ibasis, inbtypes(ibasis), re(:,ibasis,ip), ee(ibasis,ip), fe(:,ibasis,ip), &
ve(:,ibasis,ip), vele(:,ibasis,ip)
end do
else
do inod =1, n*bnum
read(11,*) ip, ibasis, inbtypes(ibasis), re(:,ibasis,ip), ee(ibasis,ip), fe(:,ibasis,ip), ve(:,ibasis,ip)
end do
end if
do i = 1, bnum
inbtypes(ibasis) = atom_type_map(inbtypes(ibasis))
end do
@ -1110,11 +1137,12 @@ module io
end if
call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re)
call add_element_data(ele_num, ee, fe, ve)
if(read_vel) vel_node(:, :, :, ele_num) = vele(:,1:max_basisnum, 1:max_ng_node)
end do
end if
call set_max_esize
return
end subroutine
end subroutine read_pycac_out
subroutine read_lmp(file, displace, temp_box_bd)
!This subroutine is used to read .cac files which are used with the lammpsCAC format
@ -1184,7 +1212,7 @@ module io
!Read atomic masses
do i = 1, type_in
read(11,*) j, mass
call ATOMMASSSPECIES(mass, atom_species)
call realatomspecies(mass, atom_species)
call add_atom_type(atom_species, type_map(i), all_new)
end do
@ -1272,7 +1300,7 @@ module io
!Read atomic masses
do i = 1, type_in
read(11,*) j, mass, textholder
call ATOMMASSSPECIES(mass, atom_species)
call realatomspecies(mass, atom_species)
call add_atom_type(atom_species, type_map(i), all_new)
end do

@ -402,7 +402,7 @@ module opt_group
!Choose what to based on what the option string is
select case(trim(textholder))
case('displace')
case('shift')
displace = .true.
do i = 1,3
arg_pos = arg_pos + 1

Loading…
Cancel
Save