module mod_crystal use mod_global use mod_lattice, only : Lattice implicit none type Crystal type(Lattice) :: LatticeType integer :: numSpecies character(len=2),allocatable :: basisAtomsSpecies(:) integer(c_int),allocatable :: basisAtomsSpeciesCode(:) real(c_double),allocatable :: basisAtomsCoords(:,:) end type Crystal contains subroutine allocateCrystal(crystalType,numBasisAtoms) implicit none type(Crystal),intent(inout) :: crystalType integer, intent(in) :: numBasisAtoms allocate(crystalType%basisAtomsSpecies(numBasisAtoms)) allocate(crystalType%basisAtomsSpeciesCode(numBasisAtoms)) allocate(crystalType%basisAtomsCoords(DIM,numBasisAtoms)) end subroutine allocateCrystal subroutine fillCrystal(crystalType,numSpecies,basisAtomsCoords,basisAtomsSpecies) implicit none type(Crystal),intent(inout) :: crystalType integer, intent(in) :: numSpecies real(c_double),intent(inout):: basisAtomsCoords(:,:) character(len=*) :: basisAtomsSpecies(:) crystalType%numSpecies = numSpecies crystalType%basisAtomsSpecies = basisAtomsSpecies crystalType%basisAtomsCoords = basisAtomsCoords end subroutine fillCrystal subroutine destroyCrystal(crystalType) implicit none type(Crystal),intent(inout) :: crystalType deallocate(crystalType%basisAtomsSpecies) deallocate(crystalType%basisAtomsCoords) deallocate(crystalType%basisAtomsSpeciesCode) end subroutine destroyCrystal subroutine constructGlobalCrystal(numCellsPerSide,cT,coords,speciesCode) implicit none integer, intent(in) :: numCellsPerSide type(Crystal), intent(in) :: cT real(c_double), intent(out) :: coords(DIM,*) integer(c_int) :: speciesCode(*) ! Local variables integer(c_int) :: i,j,k,m,a real(c_double) :: latVec(DIM),center(DIM) ! Center the crystal such the shift vectors are centered center = (numCellsPerside/2)*sum(cT%latticeType%latVec,2) + & sum(cT%basisAtomsCoords,2)/size(cT%basisAtomsCoords,2) ! make sure that crystalType is ready before proceeding ! Nikhil a = size(cT%basisAtomsCoords,2) do i=1,numCellsPerSide do j=1,numCellsPerSide do k=1,numCellsPerSide latVec = (i-1) * cT%latticeType%latVec(:,1) + & (j-1) * cT%latticeType%latVec(:,2) + & (k-1) * cT%latticeType%latVec(:,3) do m=1,size(cT%basisAtomsCoords,2) a = a+1 coords(:,a) = latVec + cT%basisAtomsCoords(:,m) - center speciesCode(a) = cT%basisAtomsSpeciesCode(m) if ((i.eq.numCellsPerside/2+1).and.(j.eq.numCellsPerSide/2+1) .and. & (k.eq.numCellsPerSide/2+1) ) then coords(:,m) = latVec + cT%basisAtomsCoords(:,m) - center! put middle particle as #1 speciesCode(m) = cT%basisAtomsSpeciesCode(m) a = a - 1 endif enddo enddo enddo enddo end subroutine constructGlobalCrystal subroutine scaleCrystal(crystalType_scaled,crystalType,scaleX,scaleY,scaleZ) implicit none ! Passed variables type(Crystal), intent(in) :: crystalType type(Crystal), intent(out) :: crystalType_scaled real(c_double), intent(in) :: scaleX,scaleY,scaleZ ! Local variables integer(c_int) :: numBasisAtoms numBasisAtoms = size(crystalType%basisAtomsCoords,2) !call allocateCrystal(scaleCrystal,numBasisAtoms) crystalType_scaled= crystalType crystalType_scaled%LatticeType%latVec(1,:) = scaleX * crystalType%LatticeType%latVec(1,:) crystalType_scaled%LatticeType%latVec(2,:) = scaleY * crystalType%LatticeType%latVec(2,:) crystalType_scaled%LatticeType%latVec(3,:) = scaleZ * crystalType%LatticeType%latVec(3,:) crystalType_scaled%basisAtomsCoords(1,:) = scaleX * crystalType%basisAtomsCoords(1,:) crystalType_scaled%basisAtomsCoords(2,:) = scaleY * crystalType%basisAtomsCoords(2,:) crystalType_scaled%basisAtomsCoords(3,:) = scaleZ * crystalType%basisAtomsCoords(3,:) crystalType_scaled%LatticeType%latVol = scaleX*scaleY*scaleZ*crystalType_scaled%LatticeType%latVol end subroutine scaleCrystal end module mod_crystal