From f0fd76f12d8d30039009507142694ceb47615dee Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 14 Sep 2021 22:33:47 -0400 Subject: [PATCH] Latest working version of cacmb --- src/atoms.f90 | 261 ++++++++++++++++++++++++++++++++++++++++++++++ src/elements.f90 | 39 ++++++- src/functions.f90 | 13 +++ src/io.f90 | 130 ++++++++++++++--------- src/opt_group.f90 | 2 +- 5 files changed, 389 insertions(+), 56 deletions(-) diff --git a/src/atoms.f90 b/src/atoms.f90 index 6da9189..dbf7956 100644 --- a/src/atoms.f90 +++ b/src/atoms.f90 @@ -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 diff --git a/src/elements.f90 b/src/elements.f90 index cb64943..6ce67b8 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -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 diff --git a/src/functions.f90 b/src/functions.f90 index ebd9a50..e0bdd16 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -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) diff --git a/src/io.f90 b/src/io.f90 index 193e10b..43a6a80 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -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 - write(11,13) + 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,21 +576,38 @@ module io case('20fcc') write(11, '(4i16)') i, basisnum(lat_ele(i)), 3, (size_ele(i)-1) end select - 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) + 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 - 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 - write(11,14) - do i = 1, atom_num - write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) - end do + 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) @@ -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(:) - read(line,*) tag, type, ra(:), ea, fa(:), va(:) + 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 - 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 + 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 diff --git a/src/opt_group.f90 b/src/opt_group.f90 index c9236da..69cb818 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -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