diff --git a/README.md b/README.md index 4714398..624d719 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,480 @@ -# VVCORElib -Projected velocity velocity autocorrelation function analyses code +# VVCORElib — GPU Branch + +**GPU-Accelerated Projected Velocity-Velocity Autocorrelation Function Analysis** + +A high-performance, GPU-accelerated Python library for computing velocity-velocity autocorrelation functions (VACF) and current correlation functions from molecular dynamics simulations. This branch leverages NVIDIA GPUs via CuPy for massively parallel computation, with MPI support for multi-GPU scaling. + +--- + +## Features + +- **GPU Acceleration** — CUDA-powered computations via CuPy for 10-100× speedups +- **Multi-GPU Support** — MPI parallelization across multiple GPUs +- **Automatic GPU Memory Management** — Smart partitioning based on available VRAM +- **Multiple Lattice Types** — Built-in support for FCC, BCC, and SC crystal structures +- **Custom Q-Grids** — Load arbitrary wavevector grids from HDF5 files +- **Current Projections** — Longitudinal (L), transverse (T), and total current correlations +- **HPC Ready** — Designed for GPU clusters with CUDA-aware MPI + +--- + +## Requirements + +### Hardware +- NVIDIA GPU with CUDA support (Compute Capability 6.0+) +- Recommended: A100, V100, or newer for optimal performance + +### Software +- Python 3.8+ +- CUDA Toolkit 11.0+ +- CUDA-aware MPI implementation (recommended) + +--- + +## Installation + +### 1. Set up a GPU-enabled environment + +```bash +# Create conda environment +conda create -n vvcore-gpu python=3.10 +conda activate vvcore-gpu + +# Install CUDA-aware mpi4py (for multi-GPU) +conda install -c conda-forge mpi4py + +# Or build from source with CUDA support: +# MPICC="mpicc" pip install mpi4py --no-cache-dir +``` + +### 2. Install dependencies + +```bash +pip install cupy-cuda11x # Match your CUDA version (cuda11x, cuda12x) +pip install numpy h5py nvidia-ml-py3 +``` + +### 3. Clone the repository + +```bash +git clone https://github.com/username/VVCORElib.git +cd VVCORElib +git checkout VVCORE_gpu +``` + +--- + +## Dependencies + +| Package | Purpose | +|----------------|--------------------------------------------| +| `cupy` | GPU array operations (NumPy API on CUDA) | +| `numpy` | CPU array operations and I/O | +| `h5py` | HDF5 file I/O for trajectories | +| `mpi4py` | MPI parallelization for multi-GPU | +| `nvidia-ml-py3`| GPU memory monitoring (`nvidia_smi`) | + +--- + +## Quick Start + +### Single GPU + +```bash +python VVCORE.py 10000 # Process 10000 frames +``` + +### Multi-GPU with MPI + +```bash +mpirun -np 4 python VVCORE.py 10000 # 4 processes across available GPUs +``` + +### HPC Cluster (SLURM + Cray MPI) + +```bash +#!/bin/bash +#SBATCH --nodes=1 +#SBATCH --gpus=4 +#SBATCH --ntasks=128 +#SBATCH --time=01:00:00 + +module load PrgEnv-nvidia cray-mpich cudatoolkit craype-accel-nvidia80 python +conda activate vvcore-gpu + +export OMP_NUM_THREADS=1 +export USE_SIMPLE_THREADED_LEVEL3=1 +export MPICH_GPU_SUPPORT_ENABLED=1 + +srun -n 128 -G 4 --cpu-bind=cores --gpu-bind=none python VVCORE.py 10000 +``` + +--- + +## Configuration + +Edit `VVCORE.py` to configure your analysis: + +```python +# Output options +opts = ['cur', 'cur_T', 'cur_L'] # Current components to compute + +# Paths +path = "./data" +traj_file = f"{path}/trajectory.h5" + +# Number of frames +N = 10000 # Or pass via command line: sys.argv[1] + +# Crystal structure +lattice = 'fcc' +a = 6.125 # Lattice constant in Å +Nq = 30 # Number of q-points along path +``` + +--- + +## Observable Options + +| Option | Description | +|----------|---------------------------------------------| +| `dens` | Density correlation function ρ(q,t) | +| `cur` | Total current j(q,t) | +| `cur_L` | Longitudinal current j_L(q,t) | +| `cur_T` | Transverse current j_T(q,t) | + +--- + +## Input Files + +### Trajectory File (HDF5) + +Trajectories must be in **H5MD** format: + +``` +trajectory.h5 +├── particles/ +│ └── all/ +│ ├── position/ +│ │ └── value [N_frames × N_atoms × 3] +│ └── velocity/ +│ └── value [N_frames × N_atoms × 3] +``` + +### Atom Type Indices (`ind.h5`) + +Required file specifying atom type indices: + +``` +ind.h5 +├── 1 [indices of type 1 atoms] +├── 2 [indices of type 2 atoms] +└── ... +``` + +**Creating an index file:** + +```python +import h5py +import numpy as np + +with h5py.File('ind.h5', 'w') as f: + f.create_dataset('1', data=np.arange(0, 500)) # Na atoms: 0-499 + f.create_dataset('2', data=np.arange(500, 1000)) # Br atoms: 500-999 +``` + +### Custom Q-Grid (`qgrid.h5`, optional) + +For arbitrary wavevector grids: + +```python +import h5py +import numpy as np + +# Custom q-points array [N_q × 3] +q_points = np.array([...]) + +with h5py.File('qgrid.h5', 'w') as f: + f.create_dataset('qx', data=q_points[:, 0]) + f.create_dataset('qy', data=q_points[:, 1]) + f.create_dataset('qz', data=q_points[:, 2]) +``` + +--- + +## Output Files + +| File | Description | +|-------------|---------------------------------------| +| `cur.h5` | Total current j(q,t) | +| `cur_L.h5` | Longitudinal current j_L(q,t) | +| `cur_T.h5` | Transverse current j_T(q,t) | +| `dens.h5` | Density ρ(q,t) (if computed) | + +Each file contains datasets keyed by atom type (e.g., `'1'`, `'2'`). + +--- + +## Python API + +### Core Modules + +| Module | Purpose | +|-----------------------|------------------------------------------------| +| `signals.py` | Main signal computation (GPU) | +| `signals_cp.py` | Alternative CuPy implementations | +| `qgrids.py` | Q-vector grid generators | +| `trajectory_reader.py`| HDF5 trajectory I/O with GPU transfer | +| `utils.py` | Memory management, MPI utilities, I/O | + +### Example: Custom Analysis + +```python +import cupy as cp +import numpy as np +from mpi4py import MPI +import nvidia_smi + +from signals import compute_signal +from utils import stack, dict_from_device, save_signal + +comm = MPI.COMM_WORLD + +# Initialize GPU +nvidia_smi.nvmlInit() +num_gpus = nvidia_smi.nvmlDeviceGetCount() +cp.cuda.Device(comm.rank % num_gpus).use() + +handle = nvidia_smi.nvmlDeviceGetHandleByIndex(comm.rank % num_gpus) +info = nvidia_smi.nvmlDeviceGetMemoryInfo(handle) + +# Compute signals +opts = ['cur', 'cur_L', 'cur_T'] +res_gpu, read_time, compute_time = compute_signal( + traj_file="./data/trajectory.h5", + N=10000, + Nq=50, + lattice='fcc', + a=4.05, + opts=opts, + comm=comm, + num_of_devices=num_gpus, + handle=handle, + info=info +) + +# Transfer from GPU to CPU +res = dict_from_device(res_gpu) +cp.cuda.get_current_stream().synchronize() + +# Gather and save +for opt in opts: + gathered = comm.gather(res[opt], root=0) + if comm.rank == 0: + stacked = stack(gathered) + save_signal(stacked, f"{opt}.h5") + +nvidia_smi.nvmlShutdown() +``` + +### Q-Grid Generators + +```python +from qgrids import grids + +# FCC high-symmetry path: Γ-X-Γ-L +Q, Qn = grids['fcc'](Nq=50, a=4.05) + +# BCC high-symmetry path: Γ-H-Γ-N +Q, Qn = grids['bcc'](Nq=50, a=2.87) + +# Simple cubic: Γ-X-M-Γ-R-M +Q, Qn = grids['sc'](Nq=50, a=3.0) + +# Load from file +Q, Qn = grids['file']() # Reads qgrid.h5 +``` + +--- + +## High-Symmetry Paths + +### FCC (Face-Centered Cubic) +``` +Γ → X → Γ → L +(0,0,0) → (0,2π/a,0) → (2π/a,2π/a,0) → (π/a,π/a,π/a) +``` + +### BCC (Body-Centered Cubic) +``` +Γ → H → Γ → N +(0,0,0) → (0,0,2π/a) → (2π/a,2π/a,2π/a) → (0,π/a,π/a) +``` + +### SC (Simple Cubic) +``` +Γ → X → M → Γ → R → M +``` + +--- + +## Theory + +### Density in Reciprocal Space + +$$\rho(\mathbf{q}, t) = \sum_{j=1}^{N} e^{i\mathbf{q} \cdot \mathbf{r}_j(t)}$$ + +### Current Density + +$$\mathbf{j}(\mathbf{q}, t) = \sum_{j=1}^{N} \mathbf{v}_j(t) \, e^{i\mathbf{q} \cdot \mathbf{r}_j(t)}$$ + +### Longitudinal Current + +$$j_L(\mathbf{q}, t) = \hat{\mathbf{q}} \cdot \mathbf{j}(\mathbf{q}, t)$$ + +### Transverse Current + +$$\mathbf{j}_T(\mathbf{q}, t) = \mathbf{j}(\mathbf{q}, t) - j_L(\mathbf{q}, t) \, \hat{\mathbf{q}}$$ + +--- + +## GPU Memory Management + +The library automatically partitions work based on available GPU memory: + +```python +# From utils.py +def signal_mem(mem, Natoms, Nq, opts): + """Returns optimal partitioning based on available VRAM""" + itemsize = 16 # complex128 + alpha = 0.50 # Use 50% of available memory + + if 'cur_T' in opts: + max_mem = 2 * itemsize * Natoms * Nq * 3 / 2**20 + else: + max_mem = 2 * itemsize * Natoms * Nq / 2**20 + + if max_mem < mem * alpha: + return int(mem * alpha) // int(max_mem), Nq + else: + return 1, int(Nq * mem * alpha / max_mem) +``` + +--- + +## Multi-GPU Scaling + +The library supports multiple MPI ranks per GPU for CPU-bound I/O overlap: + +```bash +# 4 GPUs, 32 MPI ranks (8 ranks per GPU) +mpirun -np 32 python VVCORE.py 10000 +``` + +GPU assignment is automatic: +```python +cp.cuda.Device(comm.rank % num_of_devices).use() +``` + +--- + +## Performance Tips + +1. **GPU Memory**: Larger `Nq` values require more VRAM; the library auto-partitions +2. **MPI Ranks**: Use 4-8 ranks per GPU for I/O overlap +3. **CUDA-Aware MPI**: Enable for direct GPU-GPU communication +4. **HDF5 I/O**: Keep trajectories on fast storage (NVMe, parallel FS) +5. **Batch Size**: Larger frame chunks (`Nsplit`) reduce kernel launch overhead + +--- + +## Troubleshooting + +### CUDA Out of Memory +- Reduce `Nq` or number of atoms +- Check `alpha` parameter in `signal_mem()` (default 0.5) +- Use fewer MPI ranks per GPU + +### Slow Performance +- Ensure CUDA-aware MPI is enabled +- Check GPU utilization with `nvidia-smi` +- Profile with `nvprof` or `nsys` + +### MPI Errors with GPU Arrays +- Ensure `dict_from_device()` is called before `comm.gather()` +- Synchronize streams: `cp.cuda.get_current_stream().synchronize()` + +--- + +## Comparison with MPI Branch + +| Feature | GPU Branch | MPI Branch | +|----------------------|---------------------|--------------------------| +| Compute Backend | CuPy (CUDA) | NumPy + C Extension | +| Parallelization | MPI + Multi-GPU | MPI only | +| Memory | GPU VRAM | System RAM | +| Best for | Large systems | CPU clusters | +| Autocorrelation | Not yet implemented | ✓ Implemented | +| Fourier Transform | Not yet implemented | ✓ Implemented | +| Normal Modes | Not yet implemented | ✓ Implemented | + +--- + +## Roadmap + +- [ ] GPU-accelerated autocorrelation (`compute_auto`) +- [ ] GPU-accelerated inverse Fourier transform (`compute_ift`) +- [ ] cuFFT integration for spectral analysis +- [ ] Normal mode projections on GPU +- [ ] Multi-node GPU support + +--- + +## File Structure + +``` +VVCORElib/ +├── VVCORE.py # Main entry point +├── signals.py # GPU signal computation +├── signals_cp.py # Alternative CuPy implementations +├── qgrids.py # Q-vector grid generators +├── trajectory_reader.py # HDF5 trajectory reader with GPU transfer +├── utils.py # Utilities (memory, MPI, I/O) +├── run.sh # Example SLURM/Cray submission script +├── data/ +│ └── NaBr_T300K.h5 # Example trajectory +└── ind.h5 # Atom type indices +``` + +--- + +## Citation + +If you use VVCORElib in your research, please cite: + +```bibtex +@software{vvcorelib_gpu, + title = {VVCORElib: GPU-Accelerated Velocity Autocorrelation Analysis}, + author = {Author Name}, + year = {2025}, + url = {https://github.com/username/VVCORElib}, + note = {GPU Branch} +} +``` + +--- + +## License + +[Add license information here] + +--- + +## Contributing + +Contributions are welcome! Key areas for improvement: +- GPU autocorrelation implementation +- cuFFT-based spectral analysis +- Additional lattice types (HCP) +- Performance benchmarks + diff --git a/VVCORE.py b/VVCORE.py new file mode 100644 index 0000000..833219c --- /dev/null +++ b/VVCORE.py @@ -0,0 +1,87 @@ +import time + +import cupy as cp +import numpy as np +import nvidia_smi + +from signals import compute_signal +from utils import stack, dict_from_device, save_signal +import sys + +if __name__ == "__main__": + + start_time_program = time.time() + from mpi4py import MPI + comm = MPI.COMM_WORLD + + nvidia_smi.nvmlInit() + + num_of_devices = nvidia_smi.nvmlDeviceGetCount() + cp.cuda.Device(comm.rank%num_of_devices).use() + + handle = nvidia_smi.nvmlDeviceGetHandleByIndex(comm.rank%num_of_devices) + info = nvidia_smi.nvmlDeviceGetMemoryInfo(handle) + + ############### Cur calculation ##################### + opts = ['cur', 'cur_T', 'cur_L'] + + + path = "./data" + traj_file = f"{path}/NaBr_T300K.h5" + + N = int(sys.argv[1]) + + # Initializing k points grid + + lattice = 'fcc' + a = 6.125 + Nq = 30 + + time_init = time.time() - start_time_program + if comm.rank == 0: + print(time_init) + start_time_signal = time.time() + ################# Compute signal part ############################################# + res_cp, read_time, compute_time = compute_signal(traj_file, N, Nq, lattice, a, opts, comm, num_of_devices, handle, info) + + time_signal = time.time() - start_time_signal + start_time_comm_1 = time.time() + + # Getting data from device + res = dict_from_device(res_cp) + cp.cuda.get_current_stream().synchronize() + + ################### Saving current ############################### + for opt in opts: + res_cur = comm.gather(res[opt], root = 0) + + if comm.rank == 0: + res_cur = stack(res_cur) + save_signal(res_cur, f"{path}/{opt}.h5") + del res_cur + res.pop(opt) + + + comm.Barrier() + time_comm_1 = time.time() - start_time_comm_1 + + time_program = time.time() - start_time_program + + time_proc = np.array([time_comm_1, time_signal, time_program, read_time, compute_time], dtype='d') + if comm.rank == 0: + time_total = np.zeros_like(time_proc) + else: + time_total = None + + nvidia_smi.nvmlShutdown() + + comm.Reduce([time_proc, MPI.DOUBLE], [time_total, MPI.DOUBLE], op=MPI.SUM, root=0) + + if comm.rank == 0: + time_total /= comm.size + print(f"Total time SIGNAL = {time_total[1]}") + print(f"\nTotal time READ = {time_total[-2]}") + print(f"Total time COMPUTE = {time_total[-1]}\n") + print(f"Total time PROGRAM = {time_total[2]}") + print(f"Total time COMM = {time_total[0]}") + diff --git a/VVCORE_mpi b/VVCORE_mpi deleted file mode 100644 index 7b0ab4a..0000000 --- a/VVCORE_mpi +++ /dev/null @@ -1,184 +0,0 @@ -#!/Users/temporary/anaconda3/envs/VVCORE_dynasor/bin/python -import numpy as np -import time -import h5py -import argparse - -# Custom dependences -from VVCORElib_mpi import save_signal, stack -from VVCORElib_mpi import compute_signal, compute_signal_normal_modes -from VVCORElib_mpi import compute_auto, stack_auto -from VVCORElib_mpi import compute_ift, stack_ift - -def parser(): - """Parser of the arguments""" - parser = argparse.ArgumentParser(description='VVCORE - autocorrelation function tool') - parser.add_argument('-path', type=str, default='./', help='Path to the folder with trajectory') - parser.add_argument('-i', type=str, default='data.dat', help='Input file with MD trajectory') - parser.add_argument('-opt', type=str, default='dens cur_T cur_L', help='Program option') - parser.add_argument('-N_frames', type=int, default=10000, help='Number of frames readed in one chunk') - parser.add_argument('-N_auto', type=int, default=1000, help='Number of frames to autocorrelate') - - parser.add_argument('-lattice', type=str, choices = ['hcp', 'fcc', 'bcc', 'sc', 'file'], default='fcc', help='Lattice of the crystal (for none qubic set qpoints explitely in qpoint.h5 file)') - parser.add_argument('-Nq_path', type=int, default=50, help='Number of q points to resolve (for cubic lattices) along path') - parser.add_argument('-Nq', type=int, default=50, help='Total number of q points to resolve (for cubic lattices)') - parser.add_argument('-a', type=float, default=None, help='Lattice constant for cubic lattices') - parser.add_argument('-c', type=float, default=None, help='Lattice constant for hexagonal lattices') - parser.add_argument('-M', type=str, default=None, help='Masses of the species') - - parser.add_argument('-ts', type=float, default=25, help='Timestep of the simulation in ps (10^(-12) s)') - parser.add_argument('-wmax', type=float, default=40, help='Frequency window') - parser.add_argument('-dw', type=float, default= 0.05, help='Frequency step') - - parser.add_argument('--partial', help='Option for computing all partial contributions', action='store_true') - parser.add_argument('--normal_modes', help='Option for computing projections onto normal modes', action='store_true') - - - parser.add_argument('--phase_noise', help='Option for adding noise to atomic vibrations (flips the velocity of random atoms)', action='store_true') - parser.add_argument('-noise_ind', type=int, default=0, help='Specifies to what atoms noise should be added') - parser.add_argument('-noise_ratio', type=float, default=0.0, help='Ratio of noisy atoms') - - args = parser.parse_args() - return args - -if __name__ == "__main__": - - start_time_program = time.time() - from mpi4py import MPI - comm = MPI.COMM_WORLD - - args = parser() - - path = args.path - traj_file = f"{path}/{args.i}" - opts = args.opt.split() - - N = args.N_frames - - - if not args.normal_modes: - lattice = args.lattice - Nq_path = args.Nq_path - Nq = args.Nq - a = args.a - c = args.c - - - - if not isinstance(args.M, type(None)): - M = args.M.split() - M = {str(i+1): np.float64(M[i]) for i in range(len(M))} - - partial = args.partial - - - phase_noise = args.phase_noise - noise_ind = args.noise_ind - noise_ratio = args.noise_ratio - - if noise_ind >= len(M.keys()): - raise ValueError('Index of atom is out of bounds') - - else: - opts = ['v_k'] - with h5py.File('./outfile.dispersion_relations.hdf5', 'r') as f: - Nq = f['eigenvectors_re'].shape[0]*f['eigenvectors_re'].shape[1] - - start_time_signal = time.time() - ################# Compute signal part ############################################# - if not args.normal_modes: - res, read_time, compute_time = compute_signal(traj_file, N, Nq_path, Nq, lattice, a, c, M, opts, partial, phase_noise, noise_ind, noise_ratio, comm) - else: - res, read_time, compute_time = compute_signal_normal_modes(traj_file, N, opts, comm) - - - time_signal = time.time() - start_time_signal - start_time_comm_1 = time.time() - - for opt in opts: - res_cur = comm.gather(res[opt], root = 0) - - if comm.rank == 0: - res_cur = stack(res_cur) - save_signal(res_cur, f"{path}/{opt}.h5") - del res_cur - res.pop(opt) - - - comm.Barrier() - time_comm_1 = time.time() - start_time_comm_1 - - ############## Compute auto part ############################################### - start_time_auto = time.time() - - N_auto = args.N_auto - res = compute_auto(path, N, Nq, N_auto, opts, comm) - - time_auto =time.time() - start_time_auto - - ############## Save auto ###################################### - start_time_comm_2 = time.time() - - for opt in opts: - res_cur = comm.gather(res[opt], root = 0) - if comm.rank == 0: - if comm.size > Nq: - res_cur = stack_auto(Nq, N//N_auto, res_cur) - else: - res_cur = stack(res_cur, dim = 1) - save_signal(res_cur, f"{path}/auto_{opt}.h5") - del res_cur - res.pop(opt) - - comm.Barrier() - time_comm_2 = time.time() - start_time_comm_2 - - ############## Compute ift part ############################################### - start_time_ift = time.time() - - wmax = args.wmax - dw = args.dw - ts = args.ts - - res = compute_ift(path, N_auto, Nq, ts, wmax, dw, opts, comm) - - time_ift = time.time() - start_time_ift - - - ############# Save ift ######################### - start_time_comm_3 = time.time() - - for opt in opts: - res_cur = comm.gather(res[opt], root = 0) - if comm.rank == 0: - if comm.size > Nq: - res_cur = stack_ift(Nq, int(wmax/dw), res_cur, comm) - else: - res_cur = stack(res_cur, dim = 1) - save_signal(res_cur, f"{path}/spec_{opt}.h5") - del res_cur - res.pop(opt) - - comm.Barrier() - time_comm_3 = time.time() - start_time_comm_3 - time_program = time.time() - start_time_program - - # Reduce timing - time_proc = np.array([time_comm_1 + time_comm_2 + time_comm_3, time_signal, time_auto, time_ift, time_program, read_time, compute_time], dtype='d') - - if comm.rank == 0: - time_total = np.zeros_like(time_proc) - else: - time_total = None - - comm.Reduce([time_proc, MPI.DOUBLE], [time_total, MPI.DOUBLE], op=MPI.SUM, root=0) - - if comm.rank == 0: - time_total /= comm.size - print(f"Total time SIGNAL = {time_total[1]}") - print(f"\nTotal time READ = {time_total[-2]}") - print(f"Total time COMPUTE = {time_total[-1]}\n") - print(f"Total time AUTO = {time_total[2]}") - print(f"Total time IFT = {time_total[3]}") - print(f"Total time PROGRAM = {time_total[4]}") - print(f"Total time COMM = {time_total[0]}") diff --git a/VVCORE_mpi_auto b/VVCORE_mpi_auto deleted file mode 100644 index b553600..0000000 --- a/VVCORE_mpi_auto +++ /dev/null @@ -1,62 +0,0 @@ -#!/Users/temporary/anaconda3/envs/VVCORE_dynasor/bin/python -import numpy as np -import time -import argparse - -# Custom dependences -from VVCORElib_mpi import save_signal, stack -from VVCORElib_mpi import compute_auto, stack_auto - - -def parser(): - """Parser of the arguments""" - parser = argparse.ArgumentParser(description='VVCORE - autocorrelation function tool') - parser.add_argument('-path', type=str, default='./', help='Path to the folder with trajectory') - parser.add_argument('-opt', type=str, default='dens cur_T cur_L', help='Program option') - parser.add_argument('-N_frames', type=int, default=10000, help='Number of frames readed in one chunk') - parser.add_argument('-N_auto', type=int, default=1000, help='Number of frames to autocorrelate') - - parser.add_argument('-Nq', type=int, default=50, help='Number of q points to resolve (for cubic lattices)') - args = parser.parse_args() - return args - -if __name__ == "__main__": - - start_time_program = time.time() - from mpi4py import MPI - comm = MPI.COMM_WORLD - - args = parser() - - path = args.path - opts = args.opt.split() - - N = args.N_frames - N_auto = args.N_auto - - Nq = args.Nq - - ############## Compute auto part ############################################### - start_time_auto = time.time() - - res = compute_auto(path, N, Nq, N_auto, opts, comm) - - time_auto =time.time() - start_time_auto - - ############## Save auto ###################################### - start_time_comm_2 = time.time() - - for opt in opts: - res_cur = comm.gather(res[opt], root = 0) - if comm.rank == 0: - if comm.size > Nq: - res_cur = stack_auto(Nq, N//N_auto, res_cur) - else: - res_cur = stack(res_cur, dim = 1) - save_signal(res_cur, f"{path}/auto_{opt}.h5") - del res_cur - res.pop(opt) - - comm.Barrier() - time_comm_2 = time.time() - start_time_comm_2 - time_program = time.time() - start_time_program \ No newline at end of file diff --git a/VVCORE_mpi_collect b/VVCORE_mpi_collect deleted file mode 100644 index c3a1b0b..0000000 --- a/VVCORE_mpi_collect +++ /dev/null @@ -1,56 +0,0 @@ -#!/Users/temporary/anaconda3/envs/VVCORE_dynasor/bin/python -import numpy as np -import time -import argparse - -# Custom dependences -from VVCORElib_mpi import save_signal, stack -from VVCORElib_mpi import collect_auto - - -def parser(): - """Parser of the arguments""" - parser = argparse.ArgumentParser(description='VVCORE - autocorrelation function tool') - parser.add_argument('-path', type=str, default='./', help='Path to the folder with folders') - parser.add_argument('-folders', type=str, default='1 2 3 4 5 6 7 8', help='Folders with trajectories') - parser.add_argument('-opt', type=str, default='dens cur_T cur_L', help='Program option') - - parser.add_argument('-Nq', type=int, default=50, help='Number of q points to resolve (for cubic lattices)') - args = parser.parse_args() - return args - -if __name__ == "__main__": - - start_time_program = time.time() - from mpi4py import MPI - comm = MPI.COMM_WORLD - - args = parser() - - path = args.path - folders = args.folders.split() - opts = args.opt.split() - - Nq = args.Nq - - ############## Compute auto part ############################################### - start_time_auto = time.time() - - res = collect_auto(path, folders, Nq, opts, comm) - - time_auto =time.time() - start_time_auto - - ############## Save auto ###################################### - start_time_comm_2 = time.time() - - for opt in opts: - res_cur = comm.gather(res[opt], root = 0) - if comm.rank == 0: - res_cur = stack(res_cur, dim = 1) - save_signal(res_cur, f"{path}/auto_{opt}.h5") - del res_cur - res.pop(opt) - - comm.Barrier() - time_comm_2 = time.time() - start_time_comm_2 - time_program = time.time() - start_time_program \ No newline at end of file diff --git a/VVCORE_mpi_ift b/VVCORE_mpi_ift deleted file mode 100644 index 7221512..0000000 --- a/VVCORE_mpi_ift +++ /dev/null @@ -1,68 +0,0 @@ -#!/Users/temporary/anaconda3/envs/VVCORE_dynasor/bin/python -import numpy as np -import time -import argparse - -# Custom dependences -from VVCORElib_mpi import save_signal, stack -from VVCORElib_mpi import compute_ift, stack_ift - -def parser(): - """Parser of the arguments""" - parser = argparse.ArgumentParser(description='VVCORE - autocorrelation function tool') - parser.add_argument('-path', type=str, default='./', help='Path to the folder with trajectory') - parser.add_argument('-opt', type=str, default='dens cur_T cur_L', help='Program option') - parser.add_argument('-N_auto', type=int, default=1000, help='Number of frames to autocorrelate') - - parser.add_argument('-Nq', type=int, default=50, help='Number of q points to resolve (for cubic lattices)') - - parser.add_argument('-ts', type=float, default=25, help='Timestep of the simulation in ps (10^(-12) s)') - parser.add_argument('-wmax', type=float, default=40, help='Frequency window') - parser.add_argument('-dw', type=float, default= 0.05, help='Frequency step') - args = parser.parse_args() - return args - -if __name__ == "__main__": - - start_time_program = time.time() - from mpi4py import MPI - comm = MPI.COMM_WORLD - - args = parser() - - path = args.path - opts = args.opt.split() - - N_auto = args.N_auto - - Nq = args.Nq - - ############## Compute ift part ############################################### - start_time_ift = time.time() - - wmax = args.wmax - dw = args.dw - ts = args.ts - - res = compute_ift(path, N_auto, Nq, ts, wmax, dw, opts, comm) - - time_ift = time.time() - start_time_ift - - - ############# Save ift ######################### - start_time_comm_3 = time.time() - - for opt in opts: - res_cur = comm.gather(res[opt], root = 0) - if comm.rank == 0: - if comm.size > Nq: - res_cur = stack_ift(Nq, int(wmax/dw), res_cur, comm) - else: - res_cur = stack(res_cur, dim = 1) - save_signal(res_cur, f"{path}/spec_{opt}.h5") - del res_cur - res.pop(opt) - - comm.Barrier() - time_comm_3 = time.time() - start_time_comm_3 - time_program = time.time() - start_time_program \ No newline at end of file diff --git a/VVCORE_mpi_reduce b/VVCORE_mpi_reduce deleted file mode 100644 index f6a812f..0000000 --- a/VVCORE_mpi_reduce +++ /dev/null @@ -1,89 +0,0 @@ -#!/Users/temporary/anaconda3/envs/VVCORE_dynasor/bin/python -import numpy as np -import time -import argparse - -# Custom dependences -from VVCORElib_mpi import save_signal, stack -from VVCORElib_mpi import compute_signal - -def parser(): - """Parser of the arguments""" - parser = argparse.ArgumentParser(description='VVCORE - autocorrelation function tool') - parser.add_argument('-path', type=str, default='./', help='Path to the folder with trajectory') - parser.add_argument('-i', type=str, default='data.dat', help='Input file with MD trajectory') - parser.add_argument('-opt', type=str, default='dens cur_T cur_L', help='Program option') - parser.add_argument('-N_frames', type=int, default=10000, help='Number of frames readed in one chunk') - - parser.add_argument('-lattice', type=str, choices = ['fcc', 'bcc', 'sc', 'file'], default='fcc', help='Lattice of the crystal (for none qubic set qpoints explitely in qpoint.h5 file)') - parser.add_argument('-Nq_path', type=int, default=50, help='Number of q points to resolve (for cubic lattices) along path') - parser.add_argument('-Nq', type=int, default=150, help='Total number of q points to resolve (for cubic lattices)') - parser.add_argument('-a', type=float, default=None, help='Lattice constant for cubic lattices') - parser.add_argument('-c', type=float, default=None, help='Lattice constant for hexagonal lattices') - parser.add_argument('-M', type=str, default=None, help='Masses of the species') - - parser.add_argument('--partial', help='Option for computing all partial contributions', action='store_true') - - parser.add_argument('--phase_noise', help='Option for adding noise to atomic vibrations (flips the velocity of random atoms)', action='store_true') - parser.add_argument('-noise_ind', type=int, default=0, help='Specifies to what atoms noise should be added') - parser.add_argument('-noise_ratio', type=float, default=0.0, help='Ratio of noisy atoms') - args = parser.parse_args() - return args - -if __name__ == "__main__": - - start_time_program = time.time() - from mpi4py import MPI - comm = MPI.COMM_WORLD - - args = parser() - - path = args.path - traj_file = f"{path}/{args.i}" - opts = args.opt.split() - - N = args.N_frames - - lattice = args.lattice - Nq_path = args.Nq_path - Nq = args.Nq - a = args.a - c = args.c - - if not isinstance(args.M, type(None)): - M = args.M.split() - M = {str(i+1): np.float64(M[i]) for i in range(len(M))} - - partial = args.partial - - phase_noise = args.phase_noise - noise_ind = args.noise_ind - noise_ratio = args.noise_ratio - - if noise_ind >= len(M.keys()): - raise ValueError('Index of atom is out of bounds') - - start_time_signal = time.time() - ################# Compute signal part #############################################3 - res, read_time, compute_time = compute_signal(traj_file, N, Nq_path, Nq, lattice, a, c, M, opts, partial, phase_noise, noise_ind, noise_ratio, comm) - - time_signal = time.time() - start_time_signal - - ################# Save signal part ################################ - start_time_comm_1 = time.time() - - for opt in opts: - res_cur = comm.gather(res[opt], root = 0) - - if comm.rank == 0: - res_cur = stack(res_cur) - save_signal(res_cur, f"{path}/{opt}.h5") - del res_cur - res.pop(opt) - - - comm.Barrier() - - time_comm_1 = time.time() - start_time_comm_1 - time_program = time.time() - start_time_program - diff --git a/VVCORElib_mpi/__init__.py b/VVCORElib_mpi/__init__.py deleted file mode 100644 index cd80421..0000000 --- a/VVCORElib_mpi/__init__.py +++ /dev/null @@ -1,8 +0,0 @@ -from .utils import save_signal, stack - -from .signals import compute_signal -from .signal_normal_modes import compute_signal as compute_signal_normal_modes -from .auto import compute_auto, stack_auto, collect_auto -from .ift import compute_ift, stack_ift - -from .qgrids import grids diff --git a/VVCORElib_mpi/auto.py b/VVCORElib_mpi/auto.py deleted file mode 100644 index 8af4ec3..0000000 --- a/VVCORElib_mpi/auto.py +++ /dev/null @@ -1,67 +0,0 @@ -import numpy as np -from scipy import signal -import h5py - -from .utils import get_start_end_2D, get_start_end - -def compute_auto(path, N, Nq, chunk_size, opts, comm): - assert(N % chunk_size == 0) - Nchunks = N//chunk_size - q_ind, chunk_ind = get_start_end_2D(Nq, Nchunks, comm.size, comm.rank) - res = {} - for opt in opts: - signal_f = h5py.File(f"{path}/{opt}.h5", 'r') - signal = {key: np.array(signal_f[key][chunk_ind[0]*chunk_size:(chunk_ind[-1]+1)*chunk_size, q_ind[0]:q_ind[-1]+1]) for key in signal_f.keys()} - - res[opt] = auto_transform(signal, chunk_size) - - signal_f.close() - return res - -def auto_transform(s, chunk_size): - signal_auto = {} - for t1 in s.keys(): - for t2 in s.keys(): - key1 = t1 + '_' + t2 - key2 = t2 + '_' + t1 - if key1 not in signal_auto.keys() and key2 not in signal_auto.keys(): - signal_auto[key1] = np.zeros((chunk_size, s[t1].shape[1]), dtype = np.complex128) - for k in range(s[t1].shape[0]//chunk_size): - for i in range(s[t1].shape[1]): - if len(s[t1].shape) == 3: - for m in range(s[t1].shape[2]): - signal_auto[key1][:, i] += signal.correlate(s[t1][k*chunk_size:(k+1)*chunk_size, i, m], s[t2][k*chunk_size:(k+1)*chunk_size, i, m], mode = "full")[-chunk_size:] - elif len(s[t1].shape) == 2: - signal_auto[key1][:, i] += signal.correlate(s[t1][k*chunk_size:(k+1)*chunk_size, i], s[t2][k*chunk_size:(k+1)*chunk_size, i], mode = "full")[-chunk_size:] - - return signal_auto - - -def stack_auto(Nq, Nchunks, s): - """Converts list of dicts to dict of stacked lists""" - s_final = {} - for k in s[0].keys(): - s_final[k] = np.zeros((s[0][k].shape[0], Nq), dtype = np.complex128) - for i in range(len(s)): - s_final[k][:, i%Nq] += s[i][k][:, 0] - s_final[k] /= Nchunks - return s_final - -def collect_auto(path, folders, Nq, opts, comm): - """Averages autocorrelated signal from multiple trajectories""" - q_ind = get_start_end(Nq, comm.size, comm.rank) - - res = {opt: {} for opt in opts} - - N_fold = len(folders) - for el in folders: - for opt in opts: - signal_f = h5py.File(f"{path}/{el}/auto_{opt}.h5", 'r') - if len(res[opt]) == 0: - res[opt] = {key: np.array(signal_f[key][:, q_ind[0]:q_ind[-1]+1])/(N_fold) for key in signal_f.keys()} - else: - for key in signal_f.keys(): - res[opt][key] += np.array(signal_f[key][:, q_ind[0]:q_ind[-1]+1])/(N_fold) - signal_f.close() - - return res \ No newline at end of file diff --git a/VVCORElib_mpi/ift.py b/VVCORElib_mpi/ift.py deleted file mode 100644 index 37c97c7..0000000 --- a/VVCORElib_mpi/ift.py +++ /dev/null @@ -1,60 +0,0 @@ -import numpy as np -import h5py -import psutil - -from .utils import get_num_iter, get_start_end_2D - -def ift_mem(mem, N, N_freq): - """Returns optimal partitioning on Nq and Number of frames for the process""" - itemsize = 8 # size of element in the array - alpha = 0.7 # the part of memory allowed to use - - max_mem = 4*itemsize*N*N_freq/2**20 - # Checks if there is enough memory for all Nq - if max_mem < mem*alpha: - return N_freq - else: - return int(N_freq*mem/(max_mem*alpha)) - -def compute_ift(path, N, Nq, ts, wmax, dw, opts, comm): - """Computes ift for the signal""" - freq = np.arange(0, wmax, dw) - time = np.arange(0, N)*ts/1000. - - q_ind, freq_ind = get_start_end_2D(Nq, len(freq), comm.size, comm.rank) - mem = psutil.virtual_memory().available/comm.size/2**20 - N_freq_split = ift_mem(mem, N, len(freq_ind)) - - res = {} - for opt in opts: - signal_f = h5py.File(f"{path}/auto_{opt}.h5", 'r') - signal = {key: np.array(signal_f[key][:, q_ind[0]:q_ind[-1]+1]) for key in signal_f.keys()} - - res[opt] = ift(signal, freq[freq_ind[0]:freq_ind[-1]+1], time, N_freq_split) - - signal_f.close() - return res - -def ift(s, freq, time, N_freq_split): - """Performs inverse Furior transform by simple summation""" - N_freq_iter = get_num_iter(len(freq), N_freq_split) - spec = {} - for key in s.keys(): - spec[key] = np.zeros((len(freq), s[key].shape[1]), dtype = np.complex128) - for k in range(N_freq_iter): - W, T = np.meshgrid(freq[k*N_freq_split:(k+1)*N_freq_split], time) - precomp = np.exp(-2.0j*np.pi*T*W) - for key in s.keys(): - for i in range(s[key].shape[1]): - spec[key][k*N_freq_split:(k+1)*N_freq_split, i] = np.mean(s[key][:, i][:, np.newaxis]*precomp, axis = 0) - return spec - -def stack_ift(Nq, Nfreq, s, comm): - """Converts list of dicts to dict of stacked lists in the case of 2D net""" - s_final = {} - for k in s[0].keys(): - s_final[k] = np.zeros((Nfreq, Nq), dtype = np.complex128) - for i in range(len(s)): - freq_ind = get_start_end_2D(Nq, Nfreq, comm.size, i)[1] - s_final[k][freq_ind[0]:freq_ind[-1]+1, i%Nq] = s[i][k][:, 0] - return s_final diff --git a/VVCORElib_mpi/qgrids.py b/VVCORElib_mpi/qgrids.py deleted file mode 100644 index e9897b8..0000000 --- a/VVCORElib_mpi/qgrids.py +++ /dev/null @@ -1,60 +0,0 @@ -import numpy as np -import h5py - - -def hcp(Nq, a, c): - """Generates high symmetry path for the body-centered cubic lattice""" - G_M = [[x/Nq*np.pi/a, x/Nq*np.pi/a/3**(1/2.), 0.0] for x in range(Nq+1)] - M_G = [[(Nq - x)/Nq*2*np.pi/a, 0.0, 0.0] for x in range(Nq+1)] - G_A = [[0.0, 0.0, x/Nq*np.pi/c] for x in range(Nq+1)] - A_L = [[x/Nq*np.pi/a, -x/Nq*np.pi/a/3**(1/2.), np.pi/c] for x in range(Nq+1)] - - Q = np.array(G_M + M_G + G_A + A_L, dtype = np.float64, order = 'C') + 1e-4 - Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] - return Q.transpose(), Qn - -def fcc(Nq, a): - """Generates high symmetry paths for the face-centered cubic lattice""" - G_X = [[0.0, x/Nq*2*np.pi/a, 0.0] for x in range(Nq+1)] - X_G = [[(Nq - x)/Nq*2*np.pi/a, (Nq-x)/Nq*2*np.pi/a, 0.0] for x in range(Nq+1)] - G_L = [[x/(Nq//2)*np.pi/a, x/(Nq//2)*np.pi/a, x/(Nq//2)*np.pi/a] for x in range((Nq//2)+1)] - - Q = np.array(G_X + X_G + G_L, dtype = np.float64, order = 'C') + 1e-4 - Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] - return Q.transpose(), Qn - - - -def sc(Nq, a): - """Generates high symmetry paths for the simple cubic lattice""" - G_X = [[0.0, x/Nq*np.pi/a, 0.0] for x in range(Nq+1)] - X_M = [[x/Nq*np.pi/a, np.pi/a, 0.0] for x in range(Nq+1)] - M_G = [[(Nq - x)/Nq*np.pi/a, (Nq-x)/Nq*np.pi/a, 0.0] for x in range(Nq+1)] - G_R = [[x/Nq*np.pi/a, x/Nq*np.pi/a, x/Nq*np.pi/a] for x in range(Nq+1)] - R_M = [[np.pi/a, np.pi/a, (Nq - x)/Nq*np.pi/a] for x in range(Nq+1)] - - Q = np.array(G_X + X_M + M_G + G_R + R_M, dtype = np.float64, order = 'C') + 1e-4 - - Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] - return Q.transpose(), Qn - -def bcc(Nq, a): - """Generates high symmetry path for the body-centered cubic lattice""" - G_H = [[0.0, 0.0, x/Nq*2*np.pi/a] for x in range(Nq+1)] - H_G = [[(Nq - x)/Nq*2*np.pi/a, (Nq-x)/Nq*2*np.pi/a, (Nq-x)/Nq*2*np.pi/a] for x in range(Nq+1)] - G_N = [[0.0, x/(Nq//2)*np.pi/a, x/(Nq//2)*np.pi/a] for x in range((Nq//2)+1)] - - Q = np.array(G_H + H_G + G_N, dtype = np.float64, order = 'C') + 1e-4 - Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] - return Q.transpose(), Qn - - -def from_file(): - qgrid_f = h5py.File("qgrid.h5") - Q = [qgrid_f[key] for key in qgrid_f.keys()] - Q = np.array(np.vstack(Q), order = 'C') - Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] - return Q.transpose(), Qn - - -grids = {'hcp': hcp, 'fcc': fcc, 'bcc': bcc, 'sc': sc, 'file': from_file} diff --git a/VVCORElib_mpi/signal_normal_modes.py b/VVCORElib_mpi/signal_normal_modes.py deleted file mode 100644 index 6fcd719..0000000 --- a/VVCORElib_mpi/signal_normal_modes.py +++ /dev/null @@ -1,90 +0,0 @@ -import numpy as np -import os -import psutil -import time -from .utils import get_types_ind, get_egv, gen_output_dict, get_num_iter, signal_mem, get_start_end -from .qgrids import grids -from .trajectory_reader import trajectory_h5 - -import ctypes -import distutils.sysconfig - -np_pointer = np.ctypeslib.ndpointer - -pointer_vel = np_pointer(dtype=np.float64, ndim=2, - flags='f_contiguous, aligned') -pointer_egv = np_pointer(dtype=np.complex128, ndim=2, - flags='f_contiguous, aligned') - -pointer_k = np_pointer(dtype=np.complex128, ndim=2, - flags='f_contiguous, aligned') -pointer_v_k = np_pointer(dtype=np.complex128, ndim=2, - flags='f_contiguous, aligned, writeable') - - -c_ext = ctypes.cdll.LoadLibrary(f"{os.path.dirname(__file__)}/../_rho_j_k_d{distutils.sysconfig.get_config_var('EXT_SUFFIX')}") - -c_ext.v_k.argtypes = (pointer_vel, ctypes.c_int, - pointer_k, ctypes.c_int, - pointer_v_k) - -def frame_to_signal(vel, egv): - np.require(vel, np.float64, ['F_CONTIGUOUS', 'ALIGNED']) - np.require(egv, np.float64, ['F_CONTIGUOUS', 'ALIGNED']) - - res = {'v_k': None} - Natoms = vel.shape[1] - Nq = egv.shape[1] - v_k = np.zeros((3,Nq), dtype = np.complex128, order = 'F') - c_ext.v_k(vel, Natoms, egv, Nq, v_k) - res['v_k'] = v_k.T - - return res - -def compute_signal(traj_file, N, opts, comm): - """Initializes computations for signal""" - np.random.seed(comm.rank) - - egv = get_egv() - - Nq = egv['1'].shape[1] - frame_ind = get_start_end(N, comm.size, comm.rank) - ind = get_types_ind() - traj = trajectory_h5(traj_file, ind, vel_flag = True) - mem = psutil.virtual_memory().available/comm.size/2**20 - Nsplit, Qsplit = signal_mem(mem, max([ind[key].size for key in ind]), Nq, opts) - - if Qsplit < 1: - print("WARNING: Your setup takes too much memory on a single process, job might cancel due to out of memory error") - Qsplit = 1 - - return signal(traj, egv, frame_ind, Nsplit, Qsplit, opts) - - -def signal(traj, egv, frame_ind, Nsplit, Qsplit, opts): - """Computes signal splitted according to memory requirenemnts""" - Nq = egv['1'].shape[1] - res = gen_output_dict(['1'], frame_ind.size, Nq, opts) - Q_ind = np.arange(Nq) - - Niter, Qiter = get_num_iter(frame_ind.size, Nsplit), get_num_iter(Nq, Qsplit) - Niter = frame_ind.size; Nsplit = 1 - read_time = 0 - compute_time = 0 - for i in range(frame_ind[0], frame_ind[-1]+1): - - start_time_read = time.time() - pos, vel = traj.read_one_frame(i) - - read_time += time.time() - start_time_read - - start_time_compute = time.time() - for j in range(Qiter): - cur_ind_q = Q_ind[j*Qsplit:(j+1)*Qsplit] - for k in traj.ind: - res_tmp = frame_to_signal(vel[traj.ind[k]].transpose(), egv[k][:, cur_ind_q[0]:cur_ind_q[-1]+1]) - res['v_k']['1'][i-frame_ind[0], cur_ind_q[0]:cur_ind_q[-1]+1] += res_tmp['v_k'] - compute_time += time.time() - start_time_compute - - traj.close() - return res, read_time, compute_time diff --git a/VVCORElib_mpi/signals.py b/VVCORElib_mpi/signals.py deleted file mode 100644 index a941c53..0000000 --- a/VVCORElib_mpi/signals.py +++ /dev/null @@ -1,142 +0,0 @@ -import numpy as np -import os -import psutil -import time -from .utils import get_types_ind, gen_output_dict, get_num_iter, signal_mem, get_start_end -from .qgrids import grids -from .trajectory_reader import trajectory_h5 - -import ctypes -import distutils.sysconfig - -np_pointer = np.ctypeslib.ndpointer - -pointer_pos = np_pointer(dtype=np.float64, ndim=2, - flags='f_contiguous, aligned') -pointer_vel = np_pointer(dtype=np.float64, ndim=2, - flags='f_contiguous, aligned') -pointer_k = np_pointer(dtype=np.float64, ndim=2, - flags='f_contiguous, aligned') - -pointer_dens = np_pointer(dtype=np.complex128, ndim=1, - flags='f_contiguous, aligned, writeable') -pointer_cur = np_pointer(dtype=np.complex128, ndim=2, - flags='f_contiguous, aligned, writeable') - - -c_ext = ctypes.cdll.LoadLibrary(f"{os.path.dirname(__file__)}/../_rho_j_k_d{distutils.sysconfig.get_config_var('EXT_SUFFIX')}") - -c_ext.rho_k.argtypes = (pointer_pos, ctypes.c_int, - pointer_k, ctypes.c_int, - pointer_dens) - -c_ext.rho_j_k.argtypes = (pointer_pos, pointer_vel, ctypes.c_int, - pointer_k, ctypes.c_int, - pointer_dens, pointer_cur) - -def _j_L(Qn, cur): - """Computes longtitugonal projected current""" - return (cur*Qn).sum(axis = 1)[..., np.newaxis]*Qn - -def _j_T(Qn, cur, cur_L): - """Computes transversivel projected current""" - if not isinstance(cur, type(None)) and not isinstance(cur_L, type(None)): - return cur - cur_L - elif not isinstance(cur, type(None)) and isinstance(cur_L, type(None)): - return cur - _j_L(Qn, cur) - -def frame_to_signal(pos, vel, Q, Qn, opts): - np.require(pos, np.float64, ['F_CONTIGUOUS', 'ALIGNED']) - np.require(Q, np.float64, ['F_CONTIGUOUS', 'ALIGNED']) - if not isinstance(vel, type(None)): - np.require(vel, np.float64, ['F_CONTIGUOUS', 'ALIGNED']) - - res = {'cur_L': None, 'cur': None} - Natoms = pos.shape[1] - Nq = Q.shape[1] - dens = np.zeros(Nq, dtype = np.complex128, order = 'F') - if len(opts) == 1 and opts[0] == 'dens': - c_ext.rho_k(pos, Natoms, Q, Nq, dens) - res['dens'] = dens - - if "cur" in opts or "cur_L" in opts or "cur_T" in opts: - cur = np.zeros((3,Nq), dtype = np.complex128, order = 'F') - c_ext.rho_j_k(pos, vel, Natoms, Q, Nq, dens, cur) - res['dens'] = dens - res['cur'] = cur.T - - if 'cur_L' in opts: - res['cur_L'] = _j_L(Qn, res['cur']) - if 'cur_T' in opts: - res['cur_T'] = _j_T(Qn, res['cur'], res['cur_L']) - - return res - -def compute_signal(traj_file, N, Nq_path, Nq, lattice, a, c, M, opts, partial, phase_noise, noise_ind, noise_ratio, comm): - """Initializes computations for signal""" - np.random.seed(comm.rank) - - if os.path.exists(f"qgrid.h5"): - Q, Qn = grids['file']() - else: - if lattice == 'hcp': - Q, Qn = grids[lattice](Nq_path-1, a, c) - else: - Q, Qn = grids[lattice](Nq_path-1, a) - - Q = Q[:, :Nq]; Qn = Qn[:Nq] - Nq = Q.shape[1] - frame_ind = get_start_end(N, comm.size, comm.rank) - ind = get_types_ind() - if isinstance(M, type(None)): - M = {key: np.float64(1) for key in ind} - traj = trajectory_h5(traj_file, ind, M, vel_flag = not (len(opts) == 1 and opts[0] == 'dens')) - mem = psutil.virtual_memory().available/comm.size/2**20 - Nsplit, Qsplit = signal_mem(mem, max([ind[key].size for key in ind]), Nq, opts) - - if Qsplit < 1: - print("WARNING: Your setup takes too much memory on a single process, job might cancel due to out of memory error") - Qsplit = 1 - - return signal(traj, Q, Qn, frame_ind, Nsplit, Qsplit, opts, partial, phase_noise, noise_ind, noise_ratio) - - -def signal(traj, Q, Qn, frame_ind, Nsplit, Qsplit, opts, partial, phase_noise, noise_ind, noise_ratio): - """Computes signal splitted according to memory requirenemnts""" - Nq = Q.shape[1] - if partial: - res = gen_output_dict(traj.ind.keys(), frame_ind.size, Nq, opts) - else: - res = gen_output_dict(['1'], frame_ind.size, Nq, opts) - Q_ind = np.arange(Nq) - - Niter, Qiter = get_num_iter(frame_ind.size, Nsplit), get_num_iter(Nq, Qsplit) - Niter = frame_ind.size; Nsplit = 1 - read_time = 0 - compute_time = 0 - for i in range(frame_ind[0], frame_ind[-1]+1): - - start_time_read = time.time() - pos, vel = traj.read_one_frame(i) - if phase_noise: - random_ind = np.random.randint(traj.ind[f'{noise_ind + 1}'].min(), traj.ind[f'{noise_ind + 1}'].max(), int(noise_ratio*vel.shape[0])) - vel[random_ind] *= -1 - - read_time += time.time() - start_time_read - - start_time_compute = time.time() - for j in range(Qiter): - cur_ind_q = Q_ind[j*Qsplit:(j+1)*Qsplit] - Qcur = Q[:, cur_ind_q[0]:cur_ind_q[-1]+1] - Qncur = Qn[cur_ind_q[0]:cur_ind_q[-1]+1] - for k in traj.ind: - res_tmp = frame_to_signal(pos[traj.ind[k]].transpose(), vel[traj.ind[k]].transpose(), Qcur, Qncur, opts) - for opt in opts: - if partial: - res[opt][k][i-frame_ind[0], cur_ind_q[0]:cur_ind_q[-1]+1] = traj.M[k]**(1/2.0)*res_tmp[opt] - else: - res[opt]['1'][i-frame_ind[0], cur_ind_q[0]:cur_ind_q[-1]+1] += traj.M[k]**(1/2.0)*res_tmp[opt] - compute_time += time.time() - start_time_compute - - traj.close() - return res, read_time, compute_time diff --git a/config.py b/config.py deleted file mode 100644 index 6f41825..0000000 --- a/config.py +++ /dev/null @@ -1,7 +0,0 @@ -# Let local_compiler be None in order to use the default compiler -local_compiler = None -local_linker = local_compiler -extra_compile_args = [] -extra_link_args = [] -local_link_shared = [] - diff --git a/data/NaBr_T300K.h5 b/data/NaBr_T300K.h5 new file mode 120000 index 0000000..9b22599 --- /dev/null +++ b/data/NaBr_T300K.h5 @@ -0,0 +1 @@ +/pscratch/sd/v/vladygin/NaBr_project/VVCORE/VVCORE_dipole_new/NaBr_T300K.h5 \ No newline at end of file diff --git a/device_count_test.py b/device_count_test.py new file mode 100644 index 0000000..37fe26c --- /dev/null +++ b/device_count_test.py @@ -0,0 +1,19 @@ +import cupy as cp +from mpi4py import MPI +import nvidia_smi + +comm = MPI.COMM_WORLD +print(comm.rank) + + +cp.cuda.Device(comm.rank).use() +nvidia_smi.nvmlInit() +handle = nvidia_smi.nvmlDeviceGetHandleByIndex(comm.rank) +x = cp.array([1, 2, 3]*(4-comm.rank)*10**8) +info = nvidia_smi.nvmlDeviceGetMemoryInfo(handle) +print(x.device) +print("Total memory:", info.total/2**20) +print("Free memory:", info.free/2**20) +print("Used memory:", info.used/2**20) + +nvidia_smi.nvmlShutdown() diff --git a/ind.h5 b/ind.h5 new file mode 100644 index 0000000..fcd3ade Binary files /dev/null and b/ind.h5 differ diff --git a/memory_h5py_test.py b/memory_h5py_test.py new file mode 100644 index 0000000..90d7a51 --- /dev/null +++ b/memory_h5py_test.py @@ -0,0 +1,42 @@ +import cupy as cp +import h5py +from datetime import datetime +from mpi4py import MPI +import nvidia_smi + +from trajectory_reader import trajectory_cp + +comm = MPI.COMM_WORLD +num_of_devices = cp.cuda.runtime.getDeviceCount() +if comm.rank == 0: + print(num_of_devices) +############### Cur calculation ##################### +opts = ['cur'] + +path = "/pscratch/sd/v/vladygin/VVCORE_benchmark/data/hdf5/" +traj_file = f"{path}/NaBr_T300.h5" + +N = 10000 +traj = trajectory_cp(traj_file, vel_flag = not (len(opts) == 1 and opts[0] == 'dens')) + + +cp.cuda.Device(comm.rank).use() +nvidia_smi.nvmlInit() +handle = nvidia_smi.nvmlDeviceGetHandleByIndex(comm.rank) + +if comm.rank == 0: + start_time = datetime.now() +tot = 1 +for i in range(tot): + pos, vel = traj.get_slice(0, N//num_of_devices//tot) + info = nvidia_smi.nvmlDeviceGetMemoryInfo(handle) + #if comm.rank == 0: + # print(pos.device) + # print("Total memory:", info.total/2**20) + # print("Free memory:", info.free/2**20) + # print("Used memory:", info.used/2**20) + +if comm.rank == 0: + end_time = datetime.now() + print(f"Total time of reading is {end_time - start_time}") +nvidia_smi.nvmlShutdown() diff --git a/memory_test.py b/memory_test.py new file mode 100644 index 0000000..b0b85d8 --- /dev/null +++ b/memory_test.py @@ -0,0 +1,14 @@ +import nvidia_smi + +nvidia_smi.nvmlInit() + +handle = nvidia_smi.nvmlDeviceGetHandleByIndex(0) +# card id 0 hardcoded here, there is also a call to get all available card ids, so we could iterate + +info = nvidia_smi.nvmlDeviceGetMemoryInfo(handle) + +print("Total memory:", info.total/2**20) +print("Free memory:", info.free/2**20) +print("Used memory:", info.used/2**20) + +nvidia_smi.nvmlShutdown() diff --git a/qgrids.py b/qgrids.py new file mode 100644 index 0000000..9b34ab1 --- /dev/null +++ b/qgrids.py @@ -0,0 +1,48 @@ +import numpy as np +import h5py + +def fcc(Nq, a): + """Generates high symmetry paths for the face-centered cubic lattice""" + G_X = [[0.0, x/Nq*2*np.pi/a, 0.0] for x in range(Nq+1)] + X_G = [[(Nq - x)/Nq*2*np.pi/a, (Nq-x)/Nq*2*np.pi/a, 0.0] for x in range(Nq+1)] + G_L = [[x/Nq*np.pi/a, x/Nq*np.pi/a, x/Nq*np.pi/a] for x in range(Nq+1)] + + Q = np.array(G_X + X_G + G_L, dtype = np.float64) + 1e-4 + Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] + return Q, Qn + + + +def sc(Nq, a): + """Generates high symmetry paths for the simple cubic lattice""" + G_X = [[0.0, x/Nq*np.pi/a, 0.0] for x in range(Nq+1)] + X_M = [[x/Nq*np.pi/a, np.pi/a, 0.0] for x in range(Nq+1)] + M_G = [[(Nq - x)/Nq*np.pi/a, (Nq-x)/Nq*np.pi/a, 0.0] for x in range(Nq+1)] + G_R = [[x/Nq*np.pi/a, x/Nq*pi/a, x/Nq*np.pi/a] for x in range(Nq+1)] + R_M = [[np.pi/a, np.pi/a, (Nq - x)/Nq*np.pi/a] for x in range(Nq+1)] + + Q = np.array(G_X + X_M + M_G + G_R + R_M, dtype = np.float64) + + Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] + return Q, Qn + +def bcc(Nq, a): + """Generates high symmetry path for the body-centered cubic lattice""" + G_H = [[0.0, 0.0, x/Nq*2*np.pi/a] for x in range(Nq+1)] + H_G = [[(Nq - x)/Nq*2*np.pi/a, (Nq-x)/Nq*2*np.pi/a, (Nq-x)/Nq*2*np.pi/a] for x in range(Nq+1)] + G_N = [[0.0, x/Nq*np.pi/a, x/Nq*np.pi/a] for x in range(Nq+1)] + + Q = np.array(G_H + H_G + G_N, dtype = np.float64) + 1e-4 + Qn = Q/np.linalg.norm(Q, axis = 1)[:, np.newaxis] + return Q, Qn + + +def from_file(): + qgrid_f = hdf5.File("qgrid.h5") + Q = [qgrid[key] for key in qgrid.keys()] + Q = np.vstack(Q) + Qn = Q/norm(Q, axis = 1)[:, np.newaxis] + return Q, Qn + + +grids = {'fcc': fcc, 'bcc': bcc, 'sc': sc, 'file': from_file} diff --git a/run.sh b/run.sh new file mode 100755 index 0000000..728852e --- /dev/null +++ b/run.sh @@ -0,0 +1,11 @@ +#!/bin/bash + +module load PrgEnv-nvidia cray-mpich cudatoolkit craype-accel-nvidia80 python +conda activate gpu-aware-mpi + +export OMP_NUM_THREADS=1 +export USE_SIMPLE_THREADED_LEVEL3=1 +export MPICH_GPU_SUPPORT_ENABLED=1 + + +srun -n 128 -G 4 --cpu-bind=cores --gpu-bind=none python VVCORE.py 10000 diff --git a/setup.py b/setup.py deleted file mode 100644 index 7dde51e..0000000 --- a/setup.py +++ /dev/null @@ -1,25 +0,0 @@ -#!/usr/bin/env python - -from distutils.core import setup, Extension - -rho_j_k_d_ext = Extension('_rho_j_k_d', - sources=['src/_rho_j_k.c'], - define_macros=[('RHOPREC', 'double')], - extra_compile_args=[], - extra_link_args=[], - ) - -rho_j_k_s_ext = Extension('_rho_j_k_s', - sources=['src/_rho_j_k.c'], - define_macros=[('RHOPREC', 'float')], - extra_compile_args=[], - extra_link_args=[], - ) - -setup (name = 'VVCORElib_mpi', - version = '1.0', - description = 'This is version of VVCORE with sequential reading of csv LAMMPS outputs', - packages=['VVCORElib_mpi'], - ext_modules = [rho_j_k_d_ext, rho_j_k_s_ext], - scripts=['VVCORE_mpi', 'VVCORE_mpi_reduce', 'VVCORE_mpi_auto', 'VVCORE_mpi_collect', 'VVCORE_mpi_ift'], - install_requires=['numpy', 'scipy', 'psutil', 'h5py', 'mpi4py']) diff --git a/signals.py b/signals.py new file mode 100644 index 0000000..7b0f886 --- /dev/null +++ b/signals.py @@ -0,0 +1,102 @@ +import cupy as cp +import numpy as np +import os +import time + +from utils import get_types_ind, gen_output_dict, get_num_iter, signal_mem, get_start_end +from qgrids import grids +from trajectory_reader import trajectory_h5 + +def _dens(p_proj): + """Computes projected density""" + return p_proj.sum(axis = 1) + +def _j(p_proj, vel): + """Computes full projected current""" + return (vel[:, :, cp.newaxis] * (p_proj[..., cp.newaxis])).sum(axis = 1) + +def _j_L(Qn, cur): + """Computes longtitugonal projected current""" + return (cur*Qn).sum(axis = 2)[..., cp.newaxis]*Qn + +def _j_T(Qn, cur, cur_L): + """Computes transversivel projected current""" + if not isinstance(cur, type(None)) and not isinstance(cur_L, type(None)): + return cur - cur_L + elif not isinstance(cur, type(None)) and isinstance(cur_L, type(None)): + return cur - _j_L(Qn, cur) + +def frame_to_signal(p_proj, vel, Qn, opts): + v_proj = None + + res = {'cur_L': None, 'cur': None} + if "cur" in opts or "cur_L" in opts or "cur_T" in opts: + res['cur'] = _j(p_proj, vel) + + if 'cur_L' in opts: + res['cur_L'] = _j_L(Qn, res['cur']) + if 'cur_T' in opts: + res['cur_T'] = _j_T(Qn, res['cur'], res['cur_L']) + if 'dens' in opts: + res['dens'] = _dens(p_proj) + + return res + +def compute_signal(traj_file, N, Nq, lattice, a, opts, comm, num_of_devices, handle, info): + """Initializes computations for signal""" + traj = trajectory_h5(traj_file, vel_flag = not (len(opts) == 1 and opts[0] == 'dens')) + + if os.path.exists(f"qrid.h5"): + Q, Qn = grids['file']() + else: + Q, Qn = grids[lattice](Nq, a) + + Q = cp.array(Q[:Nq]); Qn = cp.array(Qn[:Nq]) + frame_ind = get_start_end(N, comm.size, comm.rank) + ind = get_types_ind("ind.h5") + + + start_time_init = time.time() + + mem = info.free/2**20/(comm.size//num_of_devices) + + + Nsplit, Qsplit = signal_mem(mem, max([ind[key].size for key in ind]), Nq, opts) + + if Qsplit < 1: + print("WARNING: Your setup takes too much memory on a single process, job might cancel due to out of memory error") + Qsplit = 1 + res, read_time, compute_time = signal(traj, Q, Qn, ind, frame_ind, Nsplit, Qsplit, opts) + return res, read_time, compute_time + + +def signal(traj, Q, Qn, ind, frame_ind, Nsplit, Qsplit, opts): + """Computes signal splitted according to memory requirenemnts""" + Nq = Q.shape[0] + res = gen_output_dict(ind.keys(), frame_ind.size, Nq, opts) + Q_ind = np.arange(Nq) + Niter, Qiter = get_num_iter(frame_ind.size, Nsplit), get_num_iter(Nq, Qsplit) + #Niter = frame_ind.size; Nsplit = 1 + read_time = 0 + compute_time = 0 + for i in range(Niter): + cur_ind = frame_ind[Nsplit*i:Nsplit*(i+1)] + + start_time_read = time.time() + pos, vel = traj.get_slice(cur_ind[0], cur_ind[-1]+1) + read_time += time.time() - start_time_read + + start_time_compute = time.time() + for j in range(Qiter): + cur_ind_q = Q_ind[j*Qsplit:(j+1)*Qsplit] + + Qcur = Q[cur_ind_q[0]:cur_ind_q[-1]+1] + Qncur = Qn[cur_ind_q[0]:cur_ind_q[-1]+1] + p_proj = cp.exp(1.0j*(pos@Qcur.T)) + for k in ind: + res_tmp = frame_to_signal(p_proj[:, ind[k]], vel[:, ind[k]], Qncur, opts) + for opt in opts: + res[opt][k][cur_ind[0]-frame_ind[0]:cur_ind[-1]-frame_ind[0]+1, cur_ind_q[0]:cur_ind_q[-1]+1] = res_tmp[opt] + compute_time += time.time() - start_time_compute + traj.close() + return res, read_time, compute_time diff --git a/signals_cp.py b/signals_cp.py new file mode 100644 index 0000000..e294c8b --- /dev/null +++ b/signals_cp.py @@ -0,0 +1,21 @@ +import cupy as cp + +def dens(pos, Q): + """Computes projected density""" + #print(pos.shape) + #print(Q.shape) + return (cp.exp(1.0j*pos.dot(Q.T))).sum(axis = 1) + +def _j(pos, vel, Q, Qn): + """Computes full projected current""" + return cp.sum(vel[:, :, cp.newaxis] * (cp.exp(1.0j*pos.dot(Q.T))[..., cp.newaxis]), axis = 1) + +def _j_L(pos, vel, Q, Qn): + """Computes longtitugonal projected current""" + return (cp.sum(vel.dot(Qn.T)*cp.exp(1.0j*pos.dot(Q.T)), axis = 1)[..., cp.newaxis])*Qn + +def _j_T(pos, vel, Q, Qn): + """Computes transversivel projected current""" + return cp.sum((vel[:, :, cp.newaxis] - vel.dot(Qn.T)[..., cp.newaxis]*Qn) * (cp.exp(1.0j*pos.dot(Q.T))[..., cp.newaxis]), axis = 1) + +cur = {'cur': _j, 'cur_L': _j_L, 'cur_T': _j_T} diff --git a/src/_rho_j_k.c b/src/_rho_j_k.c deleted file mode 100644 index 041d342..0000000 --- a/src/_rho_j_k.c +++ /dev/null @@ -1,125 +0,0 @@ -#include - -#ifndef RHOPREC -#warning "Defaulting to double precision" -#define RHOPREC double -#endif - -void rho_k(const RHOPREC x_vec[][3], int N_x, - const RHOPREC k_vec[][3], int N_k, - RHOPREC (* restrict rho_k)[2]){ - - int x_i, k_i; - RHOPREC rho_ki_0, rho_ki_1; - register RHOPREC alpha; - - { - for(k_i=0; k_i