From 02cc28951e7f61f0e1ae18d21ebb393bf5d41948 Mon Sep 17 00:00:00 2001 From: Yuan Fengbo Date: Tue, 10 Dec 2019 15:19:25 +0800 Subject: [PATCH 1/9] replicate and disturb --- README.md | 26 +++++++++++++++++++ dpdata/system.py | 67 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) diff --git a/README.md b/README.md index 77d0dcbf8..3dcbfa335 100644 --- a/README.md +++ b/README.md @@ -152,3 +152,29 @@ Frame selection can be implemented by dpdata.LabeledSystem('OUTCAR').sub_system([0,-1]).to_deepmd_raw('dpmd_raw') ``` by which only the first and last frames are dumped to `dpmd_raw`. + +# replicate +dpdata will create a super cell of the current atom configuration. + +```python +dpdata.System('./POSCAR').replicate((1,2,3,) ) +``` +tuple(1,2,3) means don't copy atom configuration in x direction, make 2 copys in y direction, make 3 copys in z direction. + +# perturb +dpdata will disturb the box size and shape and change atom coordinates randomly. +`pert_num` means for each frame in the system ,the frames perturbed number generated. +`pert_fraction` means the realatic length the box change and the atom moves in each atom configuration +`pert_style` available options:`'uniform' 'normal' 'const'`. The distribution used to move box and atom coordinates. + + + +`uniform` will become uniform distribution on [-pert_fraction, pert_fraction) +`normal` will become normal distribution N(0,pert_fraction^2) +`const` means the box deform and relative distance atoms move in the box will be a constant pert_fraction + +The direction atoms move and box deform is random. + +```python3 +dpdata.System('./POSCAR').perturb(pert_num=3, pert_fraction=0.05, pert_style='uniform') +``` diff --git a/dpdata/system.py b/dpdata/system.py index aaba652b7..f9ae98947 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -579,6 +579,73 @@ def add_atom_names(self, atom_names): self.data['atom_names'].extend(atom_names) self.data['atom_numbs'].extend([0 for _ in atom_names]) + def replicate(self, ncopy): + tmp = System() + nframes = self.get_nframes() + data = self.data + tmp.data['atom_names'] = data['atom_names'].copy() + tmp.data['atom_numbs'] = (np.array(data['atom_numbs'].copy()) * np.prod(ncopy)).tolist() + tmp.data['atom_types'] = np.sort(np.tile(data['atom_types'].copy(),np.prod(ncopy))) + tmp.data['cells'] = data['cells'].copy() + for ii in range(3): + tmp.data['cells'][:,ii,:] *= ncopy[ii] + tmp.data['coords'] = np.tile(data['coords'].copy(),tuple(ncopy)+(1,1,1)) + + for xx in range(ncopy[0]): + for yy in range(ncopy[1]): + for zz in range(ncopy[2]): + tmp.data['coords'][xx,yy,zz,:,:,:] += xx * np.reshape(data['cells'][:,0,:], [-1,1,3])\ + + yy * np.reshape(data['cells'][:,1,:], [-1,1,3])\ + + zz * np.reshape(data['cells'][:,2,:], [-1,1,3]) + tmp.data['coords'] = np.reshape(np.transpose(tmp.data['coords'], [3,4,0,1,2,5]), (nframes, -1 , 3)) + return tmp + + + + def perturb(self, pert_num, pert_fraction, pert_style): + """ + np.random.randn return a 3*3 matrix , every element will be a sample of + standard normal distribution N(0,1) + """ + perturbed_system = System() + nframes = self.get_nframes() + for ii in range(nframes): + for jj in range(pert_num): + tmp_system = self[ii].copy() + cell_perturb_matrix = np.identity(3) + get_perturb_matrix(pert_fraction, pert_style, False) + tmp_system.data['cells'][0] = np.matmul(tmp_system.data['cells'][0],cell_perturb_matrix) + tmp_system.data['coords'][0] = np.matmul(tmp_system.data['coords'][0],cell_perturb_matrix) + for coord in tmp_system.data['coords'][0]: + atom_delta_matrix = get_perturb_matrix(pert_fraction, pert_style, True) + atom_delta_distance_vector = np.dot(np.diag(atom_delta_matrix),tmp_system.data['cells'][0]) + # print(135,atom_delta_matrix,np.linalg.norm(atom_delta_matrix),atom_delta_distance_vector,np.linalg.norm(atom_delta_distance_vector), tmp_system.data['cells'][0],) + tmp_system.data['coords'][0][ii] += atom_delta_distance_vector + perturbed_system.append(tmp_system) + return perturbed_system + +def get_perturb_matrix(pert_fraction, pert_style='uniform', diag_matrix_flag=True): + # perturb_matrix = None + delta_matrix = None + if pert_style == 'uniform': + if diag_matrix_flag is True: + delta_matrix = np.diag(np.random.rand(3) * 2 *pert_fraction - pert_fraction) + if diag_matrix_flag is False: + delta_matrix = np.tril(np.random.rand(3,3) * 2 *pert_fraction - pert_fraction) + elif pert_style == 'normal': + if diag_matrix_flag is True: + delta_matrix = np.diag(pert_fraction * np.random.randn(3)) + if diag_matrix_flag is False: + delta_matrix = np.tril(pert_fraction * np.random.randn(3,3)) + elif pert_style == 'const' : + tmp_vec = np.random.rand(3) + tmp_vec = pert_fraction * tmp_vec/np.linalg.norm(tmp_vec) + delta_matrix = np.diag(tmp_vec) + if delta_matrix is None: + RuntimeError('unsupported options pert_style={}, diag_matrix_flag={}'.format(pert_style, diag_matrix_flag)) + else: + pass + # perturb_matrix = np.identity(3) + delta_matrix + return delta_matrix class LabeledSystem (System): ''' From 4670f7c9a9f414003652b88ea0da4aee4d06b5b5 Mon Sep 17 00:00:00 2001 From: Yuan Fengbo Date: Wed, 11 Dec 2019 20:55:56 +0800 Subject: [PATCH 2/9] update doc --- README.md | 166 ++++++++++++++++++++++++-- dpdata/system.py | 85 ++++++------- tests/poscars/POSCAR.SiC | 16 +++ tests/poscars/POSCAR.SiC.replicate123 | 56 +++++++++ tests/test_replicate.py | 22 ++++ 5 files changed, 296 insertions(+), 49 deletions(-) create mode 100644 tests/poscars/POSCAR.SiC create mode 100644 tests/poscars/POSCAR.SiC.replicate123 create mode 100644 tests/test_replicate.py diff --git a/README.md b/README.md index 3dcbfa335..92196e179 100644 --- a/README.md +++ b/README.md @@ -163,18 +163,166 @@ tuple(1,2,3) means don't copy atom configuration in x direction, make 2 copys in # perturb dpdata will disturb the box size and shape and change atom coordinates randomly. -`pert_num` means for each frame in the system ,the frames perturbed number generated. -`pert_fraction` means the realatic length the box change and the atom moves in each atom configuration -`pert_style` available options:`'uniform' 'normal' 'const'`. The distribution used to move box and atom coordinates. +```python3 +dpdata.System('./POSCAR').perturb(pert_num=3, + box_pert_fraction=0.02, + atom_pert_fraction=0.01, + atom_pert_style='normal') +``` +### `pert_num` +Each frame in the input system will generate `pert_num` frames. + That means the command will return a system containing `frames of input system * pert_num` frames. + +### `box_pert_fraction` +A relative length that determines the length the box size change in each frame. It is just a fraction and doesn't have unit. Typlicaly, for cubic box with side length `a` , `a` will increase or decrease a random value from the intervel `[-a*box_pert_fraction, a*box_pert_fraction]`. -`uniform` will become uniform distribution on [-pert_fraction, pert_fraction) -`normal` will become normal distribution N(0,pert_fraction^2) -`const` means the box deform and relative distance atoms move in the box will be a constant pert_fraction +It will also change the shape of the box. That means an orthogonal box will become a non-orthogonal box after perturbing. The angle of inclination of the box is a random variable. + `box_pert_fraction` is also relating to the probability distribution function of the angle. + +See more details about how it will change the box below. -The direction atoms move and box deform is random. +### `atom_pert_fraction` +A relative length that determines the length atom moves in each frame. It is just a fraction and doesn't have unit. Typlicaly, for a cubic box with side length `a` , the mean value of the distance that atom moves is approximate `a*atom_pert_fraction`. -```python3 -dpdata.System('./POSCAR').perturb(pert_num=3, pert_fraction=0.05, pert_style='uniform') +### `atom_pert_style` +The probability distribution function used to change atom coordinates. +available options:`'uniform' 'normal' 'const'`. + +`uniform` means that if the box is a cube with side length `a`, how far atoms in the cube move is a random vector with max length `a*atom_pert_fraction ` + +`normal` means the squares of the distance atoms move are subject to a chi-square distribution (chi-square distribution can be seen as the sum of squares of normal distributed random variable. This is why we name this option as 'normal'.). +If the box is a cube with side length `a`, the mean value of the distance atom moves is `a*atom_pert_fraction ` + +`const` means that if the box is a cube with side length `a`, the distance atom moves is always `a*atom_pert_fraction `(For triclinic box, the distances are not equal.) + +The direction atoms move and box deformation is random. + +See more details about how atoms will move below. + +## The perturb details +--- +For each frame in the input system,dpdata will repeat the following steps `pert_num` times. That means the command will return a system containing `frames of input system * pert_num` frames. + +### first step: box deform, and atom moves correspondingly + +#### 1. generate a perturb matrix for box and atoms +--- +dpdata will generate a perturb matrix + +$L_1=\begin{pmatrix} +n_1+1 & n_2/2 & n_4/2 \\ +n_2/2 & n_3+1 & n_5/2 \\ +n_4/2 & n_5/2 & n_6+1 +\end{pmatrix}$ + +$n_i,i\in\{1,2,3,4,5,6\}$ are independent variables which are subject to uniform +distrubution on interval $[-box\_pert\_fraction，box\_pert\_fraction]$. +That is $n_i \sim U(-box\_pert\_fraction，box\_pert\_fraction)$ + +#### 2.box deforms +--- +The origin box matrix of this frame is defined as $B_o$, usually with lower lower triangular matrix form. That is: + +$origin\_box\_matrix \equiv B_o=\begin{pmatrix} +xx & 0 & 0 \\ +xy & yy & 0\\ +xz & yz & zz +\end{pmatrix}$ + +The box matrix after deformation will be +$deformed\_box\_matrix \equiv B_d = B_o L_1$ + +#### 3.atoms relocate in the new box. +--- +The atom coordinates in this frame will change correspondingly. +That is, +for atom with index $j, j \in\{0,1,2,3,...,atom\_num-1\}$ in the frame with coordinate $\vec r_j=(x_j, y_j, z_j)$, +the coordinate after deformation $\vec r_{jd}$ will become the matrix multiplication of $r_j$ and $L_1$. +That is: +$\vec r_{jd} = \vec r_j L_1$ + +### second step: perturb atoms randomly. + +#### 1. generate a random vector $\vec l_{j}$ for each atom +--- + +For each atom with index $j, j \in\{0,1,2,3,...,atom\_num-1\}$ in the frame with coordinate $\vec r_{jd}=(x_{jd}, y_{jd}, z_{jd})$, dpdata will generate a random vector $\vec l_{j}=(l_{jx}, l_{jy}, l_{jz})$ . + +The method used to generate $\vec l_j$ is described below. + +:tada: if pert_style is 'normal': + + +$l_{jx},l_{jy},l_{jz}$ are independent random variables which are subject to normal distribution with mean $\mu=0$ and variance $\sigma^2 = atom\_pert\_fraction^2/3$ . +That is: +$l_{jx},l_{jy},l_{jz} \sim N(0,atom\_pert\_fraction^2/3)$ + +$\Vert \vec l_{j} \Vert_2$ is the distance the atom moves, and it satisfies the following equation $\Vert \vec l_{j} \Vert_2^2=l_{jx}^2 + l_{jy}^2+l_{jz}^2$. +$3\Vert \vec l_{j} \Vert_2^2/atom\_pert\_fraction^2$ will be object to the chi-square distribution with 3 degrees of freedom. That is: +$\frac{3\Vert \vec l_{j} \Vert_2^2}{atom\_pert\_fraction^2} \sim \chi(3)$ +The expectation value of $\Vert \vec l_{j} \Vert_2$ will be `atom_pert_fraction` .That is: +$E(\Vert \vec l_{j} \Vert_2)=atom\_pert\_fraction$ + +:tada: if pert_style is 'uniform': + +$\vec l_{j}$ will be a vector point to a random point inside a 3D unit sphere and its internal space. + +That is: +$\{ \vec l_{j} \in \mathbb{R^3} | \Vert \vec l_{j} \Vert_2 \leq atom\_pert\_fraction \}$. + +The point is chosen with equal probability. That means for arbitrary two subsets of the 3D sphere and its internal space $\{(x,y,z)|x,y,z\in\mathbb{R},x^2+y^2+z^2\leq atom\_pert\_fraction^2\}$, if the 'volume' of the subsets are equal, the probability that the random point is in them are equal. + +The following method is used to generate such $\vec l_{j}$. + +random direction of equal probability +>let $\vec x$ become a 3 dimension standard normal distribution. That is +$\vec x=(x_1,x_2,x_3), x_i(i\in\{1,2,3\})\ are\ independent.$ +$x_i \sim N(0,1)$ +and +$\vec l_{j,unit\_surface}=\vec x /\Vert \vec x \Vert_2$ +Now $\vec l_{j,unit\_surface}$ point to a point located at the surface of the 3D unit sphere. + +random point of equal probability +> let $u$ become a random variable which is object to the uniform on interval [0,1]. That is: +$u \sim U(0,1)$ +and define $v$ as the 3th root of $u$ (because it is 3 dimension), That is: +$v\equiv u^{1/3}$ +Then the target vector $\vec l_{j}$ is +$\vec l_{j}=atom\_pert\_fraction \cdot v \cdot \vec l_{j,unit\_surface}$ + +:tada: if pert_style is 'const': + +$\vec l_j \equiv atom\_pert\_fraction \cdot \vec l_{j,unit\_surface}$ + +The definition of $\vec l_{j,unit\_surface}$ is described in detail above. +and it is obvious that +$\Vert \vec l_{j} \Vert_2=atom\_pert\_fraction$ +#### 2. atom moves according to $\vec l_{j}$ +--- + +After the $\vec l_j$ generated, the atom coordinates disturbed $\vec {r_{jd,disturbed}}$ will become +$\vec {r_{jd,disturbed}}=\vec r_{jd}+\vec l_{j} B_d$ + + +$B_d$ is the box matrix after deformatiion, described in first step +$\vec r_{jd}$ is the coordinate of atom j after deformatiion, described in first step. + +### the results the method System.perturb() return +--- + +At the beginning, dpdata will create an empty system. +```python +perturbed_system = System() ``` +After every single perturbed frame is generated, it will be appended to the `perturbed_system`. + +perturbed_system will contain `pert_num * frames of the input system` frames. + +$B_d$ will be used as the final box matrix in the perturbed frame(that is `System.data['cells'][0]`) +$\vec {r_{jd,disturbed}}$ will be used as the final coordinates in the perturbed frame (that is `System.data['coords'][0][j]`) +> note: 0 means the index of the frame, j means atom index, + + Finally dpdata will return `perturbed_system` as results. + diff --git a/dpdata/system.py b/dpdata/system.py index f9ae98947..bd520d848 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -583,13 +583,13 @@ def replicate(self, ncopy): tmp = System() nframes = self.get_nframes() data = self.data - tmp.data['atom_names'] = data['atom_names'].copy() - tmp.data['atom_numbs'] = (np.array(data['atom_numbs'].copy()) * np.prod(ncopy)).tolist() - tmp.data['atom_types'] = np.sort(np.tile(data['atom_types'].copy(),np.prod(ncopy))) - tmp.data['cells'] = data['cells'].copy() + tmp.data['atom_names'] = list(np.copy(data['atom_names'])) + tmp.data['atom_numbs'] = list(np.array(np.copy(data['atom_numbs'])) * np.prod(ncopy)) + tmp.data['atom_types'] = np.sort(np.tile(np.copy(data['atom_types']),np.prod(ncopy))) + tmp.data['cells'] = np.copy(data['cells']) for ii in range(3): tmp.data['cells'][:,ii,:] *= ncopy[ii] - tmp.data['coords'] = np.tile(data['coords'].copy(),tuple(ncopy)+(1,1,1)) + tmp.data['coords'] = np.tile(np.copy(data['coords']),tuple(ncopy)+(1,1,1)) for xx in range(ncopy[0]): for yy in range(ncopy[1]): @@ -600,52 +600,57 @@ def replicate(self, ncopy): tmp.data['coords'] = np.reshape(np.transpose(tmp.data['coords'], [3,4,0,1,2,5]), (nframes, -1 , 3)) return tmp - - - def perturb(self, pert_num, pert_fraction, pert_style): - """ - np.random.randn return a 3*3 matrix , every element will be a sample of - standard normal distribution N(0,1) - """ + def perturb(self, + pert_num, + box_pert_fraction, + atom_pert_fraction, + atom_pert_style='normal'): perturbed_system = System() nframes = self.get_nframes() for ii in range(nframes): for jj in range(pert_num): tmp_system = self[ii].copy() - cell_perturb_matrix = np.identity(3) + get_perturb_matrix(pert_fraction, pert_style, False) - tmp_system.data['cells'][0] = np.matmul(tmp_system.data['cells'][0],cell_perturb_matrix) - tmp_system.data['coords'][0] = np.matmul(tmp_system.data['coords'][0],cell_perturb_matrix) - for coord in tmp_system.data['coords'][0]: - atom_delta_matrix = get_perturb_matrix(pert_fraction, pert_style, True) - atom_delta_distance_vector = np.dot(np.diag(atom_delta_matrix),tmp_system.data['cells'][0]) + box_perturb_matrix = get_box_perturb_matrix(box_pert_fraction) + tmp_system.data['cells'][0] = np.matmul(tmp_system.data['cells'][0],box_perturb_matrix) + tmp_system.data['coords'][0] = np.matmul(tmp_system.data['coords'][0],box_perturb_matrix) + for kk in range(len(tmp_system.data['coords'][0])): + atom_perturb_vector = get_atom_perturb_vector(atom_pert_fraction, atom_pert_style) + atom_delta_vector = np.matmul(atom_perturb_vector,tmp_system.data['cells'][0]) # print(135,atom_delta_matrix,np.linalg.norm(atom_delta_matrix),atom_delta_distance_vector,np.linalg.norm(atom_delta_distance_vector), tmp_system.data['cells'][0],) - tmp_system.data['coords'][0][ii] += atom_delta_distance_vector + tmp_system.data['coords'][0][kk] += atom_delta_vector perturbed_system.append(tmp_system) return perturbed_system -def get_perturb_matrix(pert_fraction, pert_style='uniform', diag_matrix_flag=True): +def get_box_perturb_matrix(box_pert_fraction): + e = np.random.rand(6) * 2 *box_pert_fraction - box_pert_fraction + box_pert_matrix = np.array( + [[1+e[0], 0.5 * e[5], 0.5 * e[4]], + [0.5 * e[5], 1+e[1], 0.5 * e[3]], + [0.5 * e[4], 0.5 * e[3], 1+e[2]]] + ) + return box_pert_matrix + +def get_atom_perturb_vector(atom_pert_fraction, atom_pert_style='normal'): + """ + np.random.randn return a 3*3 matrix , every element will be a sample of + standard normal distribution N(0,1) + """ # perturb_matrix = None - delta_matrix = None - if pert_style == 'uniform': - if diag_matrix_flag is True: - delta_matrix = np.diag(np.random.rand(3) * 2 *pert_fraction - pert_fraction) - if diag_matrix_flag is False: - delta_matrix = np.tril(np.random.rand(3,3) * 2 *pert_fraction - pert_fraction) - elif pert_style == 'normal': - if diag_matrix_flag is True: - delta_matrix = np.diag(pert_fraction * np.random.randn(3)) - if diag_matrix_flag is False: - delta_matrix = np.tril(pert_fraction * np.random.randn(3,3)) - elif pert_style == 'const' : - tmp_vec = np.random.rand(3) - tmp_vec = pert_fraction * tmp_vec/np.linalg.norm(tmp_vec) - delta_matrix = np.diag(tmp_vec) - if delta_matrix is None: - RuntimeError('unsupported options pert_style={}, diag_matrix_flag={}'.format(pert_style, diag_matrix_flag)) + random_vector = None + if atom_pert_style == 'normal': + random_vector=(atom_pert_fraction/np.sqrt(3))*np.random.randn(3) + elif atom_pert_style == 'uniform': + e = np.random.randn(3) + random_unit_vector = e/np.linalg.norm(e) + v = np.random.rand(1)^(1/3) + random_vector = atom_pert_fraction*v*random_unit_vector + elif atom_pert_style == 'const' : + e = np.random.randn(3) + random_unit_vector = e/np.linalg.norm(e) + random_vector = atom_pert_fraction*random_unit_vector else: - pass - # perturb_matrix = np.identity(3) + delta_matrix - return delta_matrix + raise RuntimeError('unsupported options atom_pert_style={}'.format(atom_pert_style)) + return random_vector class LabeledSystem (System): ''' diff --git a/tests/poscars/POSCAR.SiC b/tests/poscars/POSCAR.SiC new file mode 100644 index 000000000..0722aaec2 --- /dev/null +++ b/tests/poscars/POSCAR.SiC @@ -0,0 +1,16 @@ +C4 Si4 +1.0 +4.1454462863815970e+00 0.0000000000000000e+00 0.0000000000000000e+00 +1.3000000000000000e-01 3.9594462863815970e+00 0.0000000000000000e+00 +2.5000000000000000e-01 3.7000000000000000e-01 4.0400062863815970e+00 +C Si +4 4 +cartesian + -0.0666948568 0.0604660432 2.1539831432 + 1.9907431432 -0.0864349568 0.0663397432 + 0.0885598432 2.0783031432 0.0622547432 + 2.0324731432 2.0217031432 2.0253831432 + -0.0448649568 0.1153571432 0.0540263432 + 1.9664731432 1.9313831432 -0.1405848568 + 2.0807731432 0.0065058032 1.9956031432 + 0.0632185432 1.9999431432 2.0438931432 diff --git a/tests/poscars/POSCAR.SiC.replicate123 b/tests/poscars/POSCAR.SiC.replicate123 new file mode 100644 index 000000000..ac911d463 --- /dev/null +++ b/tests/poscars/POSCAR.SiC.replicate123 @@ -0,0 +1,56 @@ +C24 Si24 +1.000000 + 4.14544629 0.00000000 0.00000000 + 0.26000000 7.91889257 0.00000000 + 0.75000000 1.11000000 12.12001886 +C Si +24 24 +Cartesian + -0.06669486 0.06046604 2.15398314 + 0.18330514 0.43046604 6.19398943 + 0.43330514 0.80046604 10.23399572 + 0.06330514 4.01991233 2.15398314 + 0.31330514 4.38991233 6.19398943 + 0.56330514 4.75991233 10.23399572 + 1.99074314 -0.08643496 0.06633974 + 2.24074314 0.28356504 4.10634603 + 2.49074314 0.65356504 8.14635232 + 2.12074314 3.87301133 0.06633974 + 2.37074314 4.24301133 4.10634603 + 2.62074314 4.61301133 8.14635232 + 0.08855984 2.07830314 0.06225474 + 0.33855984 2.44830314 4.10226103 + 0.58855984 2.81830314 8.14226732 + 0.21855984 6.03774943 0.06225474 + 0.46855984 6.40774943 4.10226103 + 0.71855984 6.77774943 8.14226732 + 2.03247314 2.02170314 2.02538314 + 2.28247314 2.39170314 6.06538943 + 2.53247314 2.76170314 10.10539572 + 2.16247314 5.98114943 2.02538314 + 2.41247314 6.35114943 6.06538943 + 2.66247314 6.72114943 10.10539572 + -0.04486496 0.11535714 0.05402634 + 0.20513504 0.48535714 4.09403263 + 0.45513504 0.85535714 8.13403892 + 0.08513504 4.07480343 0.05402634 + 0.33513504 4.44480343 4.09403263 + 0.58513504 4.81480343 8.13403892 + 1.96647314 1.93138314 -0.14058486 + 2.21647314 2.30138314 3.89942143 + 2.46647314 2.67138314 7.93942772 + 2.09647314 5.89082943 -0.14058486 + 2.34647314 6.26082943 3.89942143 + 2.59647314 6.63082943 7.93942772 + 2.08077314 0.00650580 1.99560314 + 2.33077314 0.37650580 6.03560943 + 2.58077314 0.74650580 10.07561572 + 2.21077314 3.96595209 1.99560314 + 2.46077314 4.33595209 6.03560943 + 2.71077314 4.70595209 10.07561572 + 0.06321854 1.99994314 2.04389314 + 0.31321854 2.36994314 6.08389943 + 0.56321854 2.73994314 10.12390572 + 0.19321854 5.95938943 2.04389314 + 0.44321854 6.32938943 6.08389943 + 0.69321854 6.69938943 10.12390572 diff --git a/tests/test_replicate.py b/tests/test_replicate.py new file mode 100644 index 000000000..18eba3d96 --- /dev/null +++ b/tests/test_replicate.py @@ -0,0 +1,22 @@ +import os +import numpy as np +import unittest +from context import dpdata +from comp_sys import CompSys + +class TestReplicate123(unittest.TestCase, CompSys): + def setUp (self) : + system_1_origin = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') + self.system_1 = system_1_origin.replicate((1,2,3,)) + self.system_2 = dpdata.System('poscars/POSCAR.SiC.replicate123',fmt='vasp/poscar') + self.places = 6 + +class TestReplicate123_not_change_origin(unittest.TestCase, CompSys): + def setUp (self) : + self.system_1 = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') + self.system_1.replicate((1,2,3,)) + self.system_2 = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') + self.places = 6 + +if __name__ == '__main__': + unittest.main() From ad5fa8619d74c3b6634891e16e3523581fbbc4e2 Mon Sep 17 00:00:00 2001 From: Yuan Fengbo Date: Wed, 11 Dec 2019 21:06:36 +0800 Subject: [PATCH 3/9] delete details --- README.md | 113 +++--------------------------------------------------- 1 file changed, 5 insertions(+), 108 deletions(-) diff --git a/README.md b/README.md index 92196e179..44c05e08d 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ **dpdata** is a python package for manipulating DeePMD-kit, VASP, LAMMPS data formats. dpdata only works with python 3.x. +[TOC] # Installation One can download the source code of dpdata by @@ -192,8 +193,8 @@ available options:`'uniform' 'normal' 'const'`. `uniform` means that if the box is a cube with side length `a`, how far atoms in the cube move is a random vector with max length `a*atom_pert_fraction ` -`normal` means the squares of the distance atoms move are subject to a chi-square distribution (chi-square distribution can be seen as the sum of squares of normal distributed random variable. This is why we name this option as 'normal'.). -If the box is a cube with side length `a`, the mean value of the distance atom moves is `a*atom_pert_fraction ` +`normal` means the squares of the distance atoms move are subject to a chi-square distribution with 3 degrees of freedom (chi-square distribution can be seen as the sum of squares of normal distributed random variable. This is why we name this option as 'normal'.). +If the box is a cube with side length `a`, the mean value of the distance atom moves is `a*atom_pert_fraction`. `const` means that if the box is a cube with side length `a`, the distance atom moves is always `a*atom_pert_fraction `(For triclinic box, the distances are not equal.) @@ -203,113 +204,9 @@ See more details about how atoms will move below. ## The perturb details --- -For each frame in the input system,dpdata will repeat the following steps `pert_num` times. That means the command will return a system containing `frames of input system * pert_num` frames. +See https://hackmd.io/@yeql5ephQLaGJGgFgpvIDw/rJsyux6aH for the implement of perturb -### first step: box deform, and atom moves correspondingly - -#### 1. generate a perturb matrix for box and atoms ---- -dpdata will generate a perturb matrix - -$L_1=\begin{pmatrix} -n_1+1 & n_2/2 & n_4/2 \\ -n_2/2 & n_3+1 & n_5/2 \\ -n_4/2 & n_5/2 & n_6+1 -\end{pmatrix}$ - -$n_i,i\in\{1,2,3,4,5,6\}$ are independent variables which are subject to uniform -distrubution on interval $[-box\_pert\_fraction，box\_pert\_fraction]$. -That is $n_i \sim U(-box\_pert\_fraction，box\_pert\_fraction)$ - -#### 2.box deforms ---- -The origin box matrix of this frame is defined as $B_o$, usually with lower lower triangular matrix form. That is: - -$origin\_box\_matrix \equiv B_o=\begin{pmatrix} -xx & 0 & 0 \\ -xy & yy & 0\\ -xz & yz & zz -\end{pmatrix}$ - -The box matrix after deformation will be -$deformed\_box\_matrix \equiv B_d = B_o L_1$ - -#### 3.atoms relocate in the new box. ---- -The atom coordinates in this frame will change correspondingly. -That is, -for atom with index $j, j \in\{0,1,2,3,...,atom\_num-1\}$ in the frame with coordinate $\vec r_j=(x_j, y_j, z_j)$, -the coordinate after deformation $\vec r_{jd}$ will become the matrix multiplication of $r_j$ and $L_1$. -That is: -$\vec r_{jd} = \vec r_j L_1$ - -### second step: perturb atoms randomly. - -#### 1. generate a random vector $\vec l_{j}$ for each atom ---- - -For each atom with index $j, j \in\{0,1,2,3,...,atom\_num-1\}$ in the frame with coordinate $\vec r_{jd}=(x_{jd}, y_{jd}, z_{jd})$, dpdata will generate a random vector $\vec l_{j}=(l_{jx}, l_{jy}, l_{jz})$ . - -The method used to generate $\vec l_j$ is described below. - -:tada: if pert_style is 'normal': - - -$l_{jx},l_{jy},l_{jz}$ are independent random variables which are subject to normal distribution with mean $\mu=0$ and variance $\sigma^2 = atom\_pert\_fraction^2/3$ . -That is: -$l_{jx},l_{jy},l_{jz} \sim N(0,atom\_pert\_fraction^2/3)$ - -$\Vert \vec l_{j} \Vert_2$ is the distance the atom moves, and it satisfies the following equation $\Vert \vec l_{j} \Vert_2^2=l_{jx}^2 + l_{jy}^2+l_{jz}^2$. -$3\Vert \vec l_{j} \Vert_2^2/atom\_pert\_fraction^2$ will be object to the chi-square distribution with 3 degrees of freedom. That is: -$\frac{3\Vert \vec l_{j} \Vert_2^2}{atom\_pert\_fraction^2} \sim \chi(3)$ -The expectation value of $\Vert \vec l_{j} \Vert_2$ will be `atom_pert_fraction` .That is: -$E(\Vert \vec l_{j} \Vert_2)=atom\_pert\_fraction$ - -:tada: if pert_style is 'uniform': - -$\vec l_{j}$ will be a vector point to a random point inside a 3D unit sphere and its internal space. - -That is: -$\{ \vec l_{j} \in \mathbb{R^3} | \Vert \vec l_{j} \Vert_2 \leq atom\_pert\_fraction \}$. - -The point is chosen with equal probability. That means for arbitrary two subsets of the 3D sphere and its internal space $\{(x,y,z)|x,y,z\in\mathbb{R},x^2+y^2+z^2\leq atom\_pert\_fraction^2\}$, if the 'volume' of the subsets are equal, the probability that the random point is in them are equal. - -The following method is used to generate such $\vec l_{j}$. - -random direction of equal probability ->let $\vec x$ become a 3 dimension standard normal distribution. That is -$\vec x=(x_1,x_2,x_3), x_i(i\in\{1,2,3\})\ are\ independent.$ -$x_i \sim N(0,1)$ -and -$\vec l_{j,unit\_surface}=\vec x /\Vert \vec x \Vert_2$ -Now $\vec l_{j,unit\_surface}$ point to a point located at the surface of the 3D unit sphere. - -random point of equal probability -> let $u$ become a random variable which is object to the uniform on interval [0,1]. That is: -$u \sim U(0,1)$ -and define $v$ as the 3th root of $u$ (because it is 3 dimension), That is: -$v\equiv u^{1/3}$ -Then the target vector $\vec l_{j}$ is -$\vec l_{j}=atom\_pert\_fraction \cdot v \cdot \vec l_{j,unit\_surface}$ - -:tada: if pert_style is 'const': - -$\vec l_j \equiv atom\_pert\_fraction \cdot \vec l_{j,unit\_surface}$ - -The definition of $\vec l_{j,unit\_surface}$ is described in detail above. -and it is obvious that -$\Vert \vec l_{j} \Vert_2=atom\_pert\_fraction$ -#### 2. atom moves according to $\vec l_{j}$ ---- - -After the $\vec l_j$ generated, the atom coordinates disturbed $\vec {r_{jd,disturbed}}$ will become -$\vec {r_{jd,disturbed}}=\vec r_{jd}+\vec l_{j} B_d$ - - -$B_d$ is the box matrix after deformatiion, described in first step -$\vec r_{jd}$ is the coordinate of atom j after deformatiion, described in first step. - -### the results the method System.perturb() return +## the results the method System.perturb() return --- At the beginning, dpdata will create an empty system. From 2115a0f785d85abdf513b6fc8c7029c56f83b752 Mon Sep 17 00:00:00 2001 From: Yuan Fengbo Date: Wed, 11 Dec 2019 21:08:27 +0800 Subject: [PATCH 4/9] modify doc --- README.md | 1 - 1 file changed, 1 deletion(-) diff --git a/README.md b/README.md index 44c05e08d..b0b6c8602 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,6 @@ **dpdata** is a python package for manipulating DeePMD-kit, VASP, LAMMPS data formats. dpdata only works with python 3.x. -[TOC] # Installation One can download the source code of dpdata by From 2149afc97fa5a0af64fa95fe8a5e9edb1f5de042 Mon Sep 17 00:00:00 2001 From: Yuan Fengbo Date: Fri, 13 Dec 2019 20:27:08 +0800 Subject: [PATCH 5/9] unittest for perturb && update doc --- README.md | 42 ++-------- dpdata/system.py | 86 +++++++++++++++++-- tests/poscars/POSCAR.SiC.const | 16 ++++ tests/poscars/POSCAR.SiC.normal | 16 ++++ tests/poscars/POSCAR.SiC.uniform | 16 ++++ tests/test_perturb.py | 137 +++++++++++++++++++++++++++++++ 6 files changed, 268 insertions(+), 45 deletions(-) create mode 100644 tests/poscars/POSCAR.SiC.const create mode 100644 tests/poscars/POSCAR.SiC.normal create mode 100644 tests/poscars/POSCAR.SiC.uniform create mode 100644 tests/test_perturb.py diff --git a/README.md b/README.md index b0b6c8602..d853dac14 100644 --- a/README.md +++ b/README.md @@ -165,18 +165,18 @@ tuple(1,2,3) means don't copy atom configuration in x direction, make 2 copys in dpdata will disturb the box size and shape and change atom coordinates randomly. ```python3 dpdata.System('./POSCAR').perturb(pert_num=3, - box_pert_fraction=0.02, - atom_pert_fraction=0.01, + box_pert_fraction=0.1, + atom_pert_fraction=0.05, atom_pert_style='normal') ``` ### `pert_num` Each frame in the input system will generate `pert_num` frames. - That means the command will return a system containing `frames of input system * pert_num` frames. +That means the command will return a system containing `frames of input system * pert_num` frames. ### `box_pert_fraction` -A relative length that determines the length the box size change in each frame. It is just a fraction and doesn't have unit. Typlicaly, for cubic box with side length `a` , `a` will increase or decrease a random value from the intervel `[-a*box_pert_fraction, a*box_pert_fraction]`. +A relative length that determines the length the box size change in each frame. It is just a fraction and doesn't have unit. Typlicaly, for cubic box with side length `side_length` , `side_length` will increase or decrease a random value from the intervel `[-side_length*box_pert_fraction, side_length*box_pert_fraction]`. It will also change the shape of the box. That means an orthogonal box will become a non-orthogonal box after perturbing. The angle of inclination of the box is a random variable. `box_pert_fraction` is also relating to the probability distribution function of the angle. @@ -184,41 +184,11 @@ It will also change the shape of the box. That means an orthogonal box will beco See more details about how it will change the box below. ### `atom_pert_fraction` -A relative length that determines the length atom moves in each frame. It is just a fraction and doesn't have unit. Typlicaly, for a cubic box with side length `a` , the mean value of the distance that atom moves is approximate `a*atom_pert_fraction`. +A relative length that determines the length atom moves in each frame. It is just a fraction and doesn't have unit. Typlicaly, for a cubic box with side length `side_length` , the mean value of the distance that atom moves is approximate `a*atom_pert_fraction`. ### `atom_pert_style` The probability distribution function used to change atom coordinates. -available options:`'uniform' 'normal' 'const'`. - -`uniform` means that if the box is a cube with side length `a`, how far atoms in the cube move is a random vector with max length `a*atom_pert_fraction ` - -`normal` means the squares of the distance atoms move are subject to a chi-square distribution with 3 degrees of freedom (chi-square distribution can be seen as the sum of squares of normal distributed random variable. This is why we name this option as 'normal'.). -If the box is a cube with side length `a`, the mean value of the distance atom moves is `a*atom_pert_fraction`. - -`const` means that if the box is a cube with side length `a`, the distance atom moves is always `a*atom_pert_fraction `(For triclinic box, the distances are not equal.) +available options:`'uniform' 'normal' 'const'`. The direction atoms move and box deformation is random. -See more details about how atoms will move below. - -## The perturb details ---- -See https://hackmd.io/@yeql5ephQLaGJGgFgpvIDw/rJsyux6aH for the implement of perturb - -## the results the method System.perturb() return ---- - -At the beginning, dpdata will create an empty system. -```python -perturbed_system = System() -``` -After every single perturbed frame is generated, it will be appended to the `perturbed_system`. - -perturbed_system will contain `pert_num * frames of the input system` frames. - -$B_d$ will be used as the final box matrix in the perturbed frame(that is `System.data['cells'][0]`) -$\vec {r_{jd,disturbed}}$ will be used as the final coordinates in the perturbed frame (that is `System.data['coords'][0][j]`) -> note: 0 means the index of the frame, j means atom index, - - Finally dpdata will return `perturbed_system` as results. - diff --git a/dpdata/system.py b/dpdata/system.py index bd520d848..e5c15e13c 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -580,6 +580,30 @@ def add_atom_names(self, atom_names): self.data['atom_numbs'].extend([0 for _ in atom_names]) def replicate(self, ncopy): + """ + Replicate the each frame in the system in 3 dimensions. + Each frame in the system will become a supercell. + + Parameters + ---------- + ncopy : + list: [4,2,3] + or tuple: (4,2,3,) + make `ncopy[0]` copys in x dimensions, + make `ncopy[1]` copys in y dimensions, + make `ncopy[2]` copys in z dimensions. + + Returns + ------- + tmp : System + The system after replication. + """ + if len(ncopy) !=3: + raise RuntimeError('ncopy must be a list or tuple with 3 int') + for ii in ncopy: + if type(ncopy(ii)) is not int: + raise RuntimeError('ncopy must be a list or tuple must with 3 int') + tmp = System() nframes = self.get_nframes() data = self.data @@ -605,6 +629,43 @@ def perturb(self, box_pert_fraction, atom_pert_fraction, atom_pert_style='normal'): + """ + Perturb each frame in the system randomly. + The box will be changed randomly, + and atoms will move random distance in random direction. + + Parameters + ---------- + pert_num : int + Each frame in the system will make `pert_num` copies, + and all the copies will be perturbed. + That means the system to be returned will contain `pert_num` * frame_num of the input system. + box_pert_fraction : float + A fraction determines how much will box deform, typically less than 0.3. + For a cubic box with side length `side_length`, + `side_length` will increase or decrease with the max length `box_pert_fraction*side_length` . + If box_pert_fraction is not zero, the shape of the box will also be changed, + and that means a orthogonal box will become a non-orthogonal box. + atom_pert_fraction : float + A fraction determines how far will atoms move, typically less than 0.3. + For a cubic box with side length `side_length`, + atoms will move about `atom_pert_fraction*side_length` in random direction. + The distribution of the distance atoms move is also determined by atom_pert_style + atom_pert_fraction : str + Determines the distribution of the distance atoms move is subject to. + Avaliable options are + - `'normal'`: the `distance` will be object to `chi-square distribution with 3 degrees of freedom` after normalization. + The mean value of the distance is about `atom_pert_fraction*side_length` + - `'uniform'`: will generate uniformly random points in a 3D-balls and transform these points linearly as vectors to be used by atoms. + And the max length of the distance atoms move is about `atom_pert_fraction*side_length` + - `'const'`: The distance atoms move will be a constant `atom_pert_fraction*side_length` for cubic box. + And for other shape boxes, the distance will not be constant + + Returns + ------- + perturbed_system : System + The perturbed_system. It contains `pert_num` * frame_num of the input system frames. + """ perturbed_system = System() nframes = self.get_nframes() for ii in range(nframes): @@ -616,13 +677,16 @@ def perturb(self, for kk in range(len(tmp_system.data['coords'][0])): atom_perturb_vector = get_atom_perturb_vector(atom_pert_fraction, atom_pert_style) atom_delta_vector = np.matmul(atom_perturb_vector,tmp_system.data['cells'][0]) - # print(135,atom_delta_matrix,np.linalg.norm(atom_delta_matrix),atom_delta_distance_vector,np.linalg.norm(atom_delta_distance_vector), tmp_system.data['cells'][0],) tmp_system.data['coords'][0][kk] += atom_delta_vector + tmp_system.rot_lower_triangular() perturbed_system.append(tmp_system) return perturbed_system def get_box_perturb_matrix(box_pert_fraction): - e = np.random.rand(6) * 2 *box_pert_fraction - box_pert_fraction + if box_pert_fraction<0: + raise RuntimeError('box_pert_fraction can not be negative') + e0 = np.random.rand(6) + e = e0 * 2 *box_pert_fraction - box_pert_fraction box_pert_matrix = np.array( [[1+e[0], 0.5 * e[5], 0.5 * e[4]], [0.5 * e[5], 1+e[1], 0.5 * e[3]], @@ -631,21 +695,25 @@ def get_box_perturb_matrix(box_pert_fraction): return box_pert_matrix def get_atom_perturb_vector(atom_pert_fraction, atom_pert_style='normal'): - """ - np.random.randn return a 3*3 matrix , every element will be a sample of - standard normal distribution N(0,1) - """ - # perturb_matrix = None random_vector = None + if atom_pert_fraction < 0: + raise RuntimeError('atom_pert_fraction can not be negative') + if atom_pert_style == 'normal': - random_vector=(atom_pert_fraction/np.sqrt(3))*np.random.randn(3) + e = np.random.randn(3) + random_vector=(atom_pert_fraction/np.sqrt(3))*e elif atom_pert_style == 'uniform': e = np.random.randn(3) + while np.linalg.norm(e) < 0.1: + e = np.random.randn(3) random_unit_vector = e/np.linalg.norm(e) - v = np.random.rand(1)^(1/3) + v0 = np.random.rand(1) + v = np.power(v0,1/3) random_vector = atom_pert_fraction*v*random_unit_vector elif atom_pert_style == 'const' : e = np.random.randn(3) + while np.linalg.norm(e) < 0.1: + e = np.random.randn(3) random_unit_vector = e/np.linalg.norm(e) random_vector = atom_pert_fraction*random_unit_vector else: diff --git a/tests/poscars/POSCAR.SiC.const b/tests/poscars/POSCAR.SiC.const new file mode 100644 index 000000000..8192f64ee --- /dev/null +++ b/tests/poscars/POSCAR.SiC.const @@ -0,0 +1,16 @@ +C4 Si4 +1.0 +3.9821590818265662e+00 -6.3518285952025531e-18 -1.3636822777448101e-18 +1.4441625743554549e-02 3.8797851714535434e+00 -3.8802125261354576e-18 +2.7058404544835651e-01 5.0846848186579174e-01 4.0806630199703786e+00 +C Si +4 4 +Cartesian + -0.1284024102 0.0411632429 2.5392886616 + 1.6605207719 -0.0253819250 0.3618554716 + 0.1640583101 2.1563760674 0.4480081454 + 2.1989105171 1.7886495520 2.0950553386 + -0.1078475228 0.2108036732 -0.3228714337 + 1.8447569702 2.1025290971 -0.4494344734 + 1.6511638763 -0.0784401477 1.9721672179 + -0.0382262924 1.6485487817 1.9674381019 diff --git a/tests/poscars/POSCAR.SiC.normal b/tests/poscars/POSCAR.SiC.normal new file mode 100644 index 000000000..cb047ce5f --- /dev/null +++ b/tests/poscars/POSCAR.SiC.normal @@ -0,0 +1,16 @@ +C4 Si4 +1.0 +4.0972699470293215e+00 1.0149069314692106e-18 -1.2199095357777910e-17 +1.0494335660165277e-01 4.0313628735769402e+00 -4.8506301604322156e-19 +4.3100124938951284e-01 4.1459063237657107e-01 4.1291342816555661e+00 +C Si +4 4 +Cartesian + -0.3125685718 0.3578306188 2.6198865483 + 2.0773182946 -0.1374738727 0.1909202951 + -0.0358416028 2.2809266017 0.5372668691 + 2.1005017067 2.2434342103 2.4385806771 + -0.1565383829 0.1676212247 -0.0494987290 + 1.7296257122 1.9113469908 -0.0379209714 + 2.4121509301 0.3178978635 2.4128037669 + -0.0098316038 1.8287827235 2.0661147519 diff --git a/tests/poscars/POSCAR.SiC.uniform b/tests/poscars/POSCAR.SiC.uniform new file mode 100644 index 000000000..7971b3653 --- /dev/null +++ b/tests/poscars/POSCAR.SiC.uniform @@ -0,0 +1,16 @@ +C4 Si4 +1.0 +4.0393044668683302e+00 3.7866840715513190e-18 -3.4832403817178164e-18 +2.4199759881076477e-01 4.1312396844887331e+00 2.2713212548364351e-17 +4.4481025327383367e-01 2.0732581915669496e-01 4.0512745001512140e+00 +C Si +4 4 +Cartesian + -0.2187320173 -0.3278268312 1.9808168543 + 1.7325893560 -0.2907437704 -0.1316407092 + 0.1697204408 2.0456083784 -0.0881243869 + 2.2762667222 2.3246948417 1.8605727611 + -0.0627770396 0.1611339796 0.0659268519 + 2.0190368430 1.8199315498 -0.3333978387 + 2.2452379290 -0.0729729887 2.3597409002 + 0.0376749365 1.7102483986 2.1298706088 diff --git a/tests/test_perturb.py b/tests/test_perturb.py new file mode 100644 index 000000000..50ca295fa --- /dev/null +++ b/tests/test_perturb.py @@ -0,0 +1,137 @@ +import os +import numpy as np +import unittest +from context import dpdata +from comp_sys import CompSys + +from unittest.mock import Mock +from unittest.mock import patch, MagicMock + +class NormalGenerator(object): + def __init__(self): + self.randn_generator = self.get_randn_generator() + self.rand_generator = self.get_rand_generator() + def randn(self,number): + return next(self.randn_generator) + def rand(self,number): + return next(self.rand_generator) + @staticmethod + def get_randn_generator(): + data = np.asarray([ + [-1.67097394, 1.00565478, 1.75499669], + [ 0.40123822, -0.26832357, 0.51643976], + [-0.70325165, 0.50154408, 1.98677553], + [-0.13296745, 0.55434201, 1.54581177], + [-0.44232103, 0.25854521, -0.43925698], + [-0.86827367, -0.27676856, 0.44365476], + [0.92187855, 1.09601508, 1.56535943], + [-0.61866517, -0.96390108, -0.09593058]]) + count = 0 + while True: + yield data[count] + count +=1 + + @staticmethod + def get_rand_generator(): + yield np.asarray([0.38114487,0.68166018,0.72880982,0.59149848,0.9529526,0.4403747 ]) + +class UniformGenerator(object): + def __init__(self): + self.randn_generator = self.get_randn_generator() + self.rand_generator = self.get_rand_generator() + def randn(self,number): + return next(self.randn_generator) + def rand(self,number): + return next(self.rand_generator) + + @staticmethod + def get_randn_generator(): + data = [[-0.6338632, -0.80482399, -0.51334222], + [-1.18911393, -1.24521526, -1.34214803], + [ 0.20640311, -0.52666381, -0.71929358], + [ 0.51478689, 1.11296436, -0.61690281], + [-0.77160527, 1.11681373, 0.31440837], + [ 0.71732744, -1.56937797, -1.60701756], + [ 0.6121493, -0.06944649, 2.58409679], + [-0.86858833, -1.44664272, 0.40813879]] + yield np.asarray([0.0001,0.0001,0.0001]) # test for not using small vector + count = 0 + while True: + yield data[count] + count +=1 + + @staticmethod + def get_rand_generator(): + data = np.asarray([[0.96216328], [0.50364801], + [0.10566744], [0.82063694], + [0.00212345], [0.3384801], + [0.75329772], [0.59874454]]) + + yield np.asarray([0.24009826,0.9346383,0.54357125,0.05932848,0.96751515,0.78645846]) + count =0 + while True: + yield np.asarray(data[count]) + count+=1 + +class ConstGenerator(object): + def __init__(self): + self.randn_generator = self.get_randn_generator() + self.rand_generator = self.get_rand_generator() + def randn(self,number): + return next(self.randn_generator) + def rand(self,number): + return next(self.rand_generator) + + @staticmethod + def get_randn_generator(): + data = [[-0.37656207, -0.50545915, 1.23769463], + [-0.77863147, 0.05874498, 0.81611906], + [0.557889, 0.36082455, 1.89695245], + [ 1.07846313, -1.07323247, 0.18513304], + [-0.19418326, 0.75757126, -1.90317556], + [ 0.16438486, 1.40294282, -1.6183854 ], + [-1.71033533, -0.73018772, -0.19946578], + [-0.29332031, -2.00584426, -0.49622181]] + yield np.asarray([0.0001,0.0001,0.0001]) # test for not using small vector + count = 0 + while True: + yield data[count] + count +=1 + + @staticmethod + def get_rand_generator(): + yield np.asarray([0.10497343,0.30020419,0.60592304, 0.85344648, 0.60080211, 0.22293347]) + +# %% +class TestPerturbNormal(unittest.TestCase, CompSys): + @patch('numpy.random') + def setUp (self, random_mock): + random_mock.rand = NormalGenerator().rand + random_mock.randn = NormalGenerator().randn + system_1_origin = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') + self.system_1 = system_1_origin.perturb(1,0.05,0.1,'normal') + self.system_2 = dpdata.System('poscars/POSCAR.SiC.normal',fmt='vasp/poscar') + self.places = 6 + +class TestPerturbUniform(unittest.TestCase, CompSys): + @patch('numpy.random') + def setUp (self, random_mock) : + random_mock.rand = UniformGenerator().rand + random_mock.randn = UniformGenerator().randn + system_1_origin = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') + self.system_1 = system_1_origin.perturb(1,0.05,0.1,'uniform') + self.system_2 = dpdata.System('poscars/POSCAR.SiC.uniform',fmt='vasp/poscar') + self.places = 6 + +class TestPerturbConst(unittest.TestCase, CompSys): + @patch('numpy.random') + def setUp (self, random_mock) : + random_mock.rand = ConstGenerator().rand + random_mock.randn = ConstGenerator().randn + system_1_origin = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') + self.system_1 = system_1_origin.perturb(1,0.05,0.1,'const') + self.system_2 = dpdata.System('poscars/POSCAR.SiC.const',fmt='vasp/poscar') + self.places = 6 + +if __name__ == '__main__': + unittest.main() \ No newline at end of file From 091948e3041bb44f3ce8c44a9505796bda9d5392 Mon Sep 17 00:00:00 2001 From: Yuan Fengbo Date: Fri, 13 Dec 2019 20:31:16 +0800 Subject: [PATCH 6/9] update replicate --- dpdata/system.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpdata/system.py b/dpdata/system.py index e5c15e13c..c02e90d4f 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -601,7 +601,7 @@ def replicate(self, ncopy): if len(ncopy) !=3: raise RuntimeError('ncopy must be a list or tuple with 3 int') for ii in ncopy: - if type(ncopy(ii)) is not int: + if type(ii) is not int: raise RuntimeError('ncopy must be a list or tuple must with 3 int') tmp = System() From ad0e99a82735239f1b72d862576cc95de4ca609d Mon Sep 17 00:00:00 2001 From: Yuan Fengbo Date: Fri, 13 Dec 2019 21:51:56 +0800 Subject: [PATCH 7/9] change relative to abs distances for atoms perturbation --- README.md | 22 +++++------ dpdata/system.py | 44 ++++++++++----------- tests/poscars/POSCAR.SiC.const | 22 +++++------ tests/poscars/POSCAR.SiC.normal | 22 +++++------ tests/poscars/POSCAR.SiC.uniform | 22 +++++------ tests/test_perturb.py | 68 ++++++++++++++++---------------- 6 files changed, 98 insertions(+), 102 deletions(-) diff --git a/README.md b/README.md index d853dac14..85f1f3387 100644 --- a/README.md +++ b/README.md @@ -163,11 +163,12 @@ tuple(1,2,3) means don't copy atom configuration in x direction, make 2 copys in # perturb dpdata will disturb the box size and shape and change atom coordinates randomly. -```python3 -dpdata.System('./POSCAR').perturb(pert_num=3, - box_pert_fraction=0.1, - atom_pert_fraction=0.05, +```python +perturbed_system = dpdata.System('./POSCAR').perturb(pert_num=3, + box_pert_fraction=0.05, + atom_pert_distance=0.6, atom_pert_style='normal') +print(perturbed_system.data) ``` ### `pert_num` @@ -179,16 +180,13 @@ That means the command will return a system containing `frames of input system * A relative length that determines the length the box size change in each frame. It is just a fraction and doesn't have unit. Typlicaly, for cubic box with side length `side_length` , `side_length` will increase or decrease a random value from the intervel `[-side_length*box_pert_fraction, side_length*box_pert_fraction]`. It will also change the shape of the box. That means an orthogonal box will become a non-orthogonal box after perturbing. The angle of inclination of the box is a random variable. - `box_pert_fraction` is also relating to the probability distribution function of the angle. - -See more details about how it will change the box below. +`box_pert_fraction` is also relating to the probability distribution function of the angle. -### `atom_pert_fraction` -A relative length that determines the length atom moves in each frame. It is just a fraction and doesn't have unit. Typlicaly, for a cubic box with side length `side_length` , the mean value of the distance that atom moves is approximate `a*atom_pert_fraction`. +### `atom_pert_distance` +unit:Angstrom. Determine the distance atoms move in each frame. +The mean value of the distance that atom moves is about `atom_pert_distance`. ### `atom_pert_style` The probability distribution function used to change atom coordinates. -available options:`'uniform' 'normal' 'const'`. - +available options:`'uniform', 'normal', 'const'`. The direction atoms move and box deformation is random. - diff --git a/dpdata/system.py b/dpdata/system.py index c02e90d4f..ab020abc5 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -627,7 +627,7 @@ def replicate(self, ncopy): def perturb(self, pert_num, box_pert_fraction, - atom_pert_fraction, + atom_pert_distance, atom_pert_style='normal'): """ Perturb each frame in the system randomly. @@ -643,23 +643,22 @@ def perturb(self, box_pert_fraction : float A fraction determines how much will box deform, typically less than 0.3. For a cubic box with side length `side_length`, - `side_length` will increase or decrease with the max length `box_pert_fraction*side_length` . + `side_length` will increase or decrease with the max value `box_pert_fraction*side_length` . If box_pert_fraction is not zero, the shape of the box will also be changed, - and that means a orthogonal box will become a non-orthogonal box. - atom_pert_fraction : float - A fraction determines how far will atoms move, typically less than 0.3. - For a cubic box with side length `side_length`, - atoms will move about `atom_pert_fraction*side_length` in random direction. - The distribution of the distance atoms move is also determined by atom_pert_style - atom_pert_fraction : str + and that means a orthogonal box will become a non-orthogonal one. + atom_pert_distance: float + unit: Angstrom. A distance determines how far atoms will move. + Atoms will move about `atom_pert_distance` in random direction. + The distribution of the distance atoms move is determined by atom_pert_style + atom_pert_style : str Determines the distribution of the distance atoms move is subject to. Avaliable options are - `'normal'`: the `distance` will be object to `chi-square distribution with 3 degrees of freedom` after normalization. - The mean value of the distance is about `atom_pert_fraction*side_length` - - `'uniform'`: will generate uniformly random points in a 3D-balls and transform these points linearly as vectors to be used by atoms. - And the max length of the distance atoms move is about `atom_pert_fraction*side_length` - - `'const'`: The distance atoms move will be a constant `atom_pert_fraction*side_length` for cubic box. - And for other shape boxes, the distance will not be constant + The mean value of the distance is `atom_pert_fraction*side_length` + - `'uniform'`: will generate uniformly random points in a 3D-balls with radius as `atom_pert_distance`. + These points are treated as vector used by atoms to move. + Obviously, the max length of the distance atoms move is `atom_pert_distance`. + - `'const'`: The distance atoms move will be a constant `atom_pert_distance`. Returns ------- @@ -675,9 +674,8 @@ def perturb(self, tmp_system.data['cells'][0] = np.matmul(tmp_system.data['cells'][0],box_perturb_matrix) tmp_system.data['coords'][0] = np.matmul(tmp_system.data['coords'][0],box_perturb_matrix) for kk in range(len(tmp_system.data['coords'][0])): - atom_perturb_vector = get_atom_perturb_vector(atom_pert_fraction, atom_pert_style) - atom_delta_vector = np.matmul(atom_perturb_vector,tmp_system.data['cells'][0]) - tmp_system.data['coords'][0][kk] += atom_delta_vector + atom_perturb_vector = get_atom_perturb_vector(atom_pert_distance, atom_pert_style) + tmp_system.data['coords'][0][kk] += atom_perturb_vector tmp_system.rot_lower_triangular() perturbed_system.append(tmp_system) return perturbed_system @@ -694,14 +692,14 @@ def get_box_perturb_matrix(box_pert_fraction): ) return box_pert_matrix -def get_atom_perturb_vector(atom_pert_fraction, atom_pert_style='normal'): +def get_atom_perturb_vector(atom_pert_distance, atom_pert_style='normal'): random_vector = None - if atom_pert_fraction < 0: - raise RuntimeError('atom_pert_fraction can not be negative') + if atom_pert_distance < 0: + raise RuntimeError('atom_pert_distance can not be negative') if atom_pert_style == 'normal': e = np.random.randn(3) - random_vector=(atom_pert_fraction/np.sqrt(3))*e + random_vector=(atom_pert_distance/np.sqrt(3))*e elif atom_pert_style == 'uniform': e = np.random.randn(3) while np.linalg.norm(e) < 0.1: @@ -709,13 +707,13 @@ def get_atom_perturb_vector(atom_pert_fraction, atom_pert_style='normal'): random_unit_vector = e/np.linalg.norm(e) v0 = np.random.rand(1) v = np.power(v0,1/3) - random_vector = atom_pert_fraction*v*random_unit_vector + random_vector = atom_pert_distance*v*random_unit_vector elif atom_pert_style == 'const' : e = np.random.randn(3) while np.linalg.norm(e) < 0.1: e = np.random.randn(3) random_unit_vector = e/np.linalg.norm(e) - random_vector = atom_pert_fraction*random_unit_vector + random_vector = atom_pert_distance*random_unit_vector else: raise RuntimeError('unsupported options atom_pert_style={}'.format(atom_pert_style)) return random_vector diff --git a/tests/poscars/POSCAR.SiC.const b/tests/poscars/POSCAR.SiC.const index 8192f64ee..986356eb4 100644 --- a/tests/poscars/POSCAR.SiC.const +++ b/tests/poscars/POSCAR.SiC.const @@ -1,16 +1,16 @@ C4 Si4 1.0 -3.9821590818265662e+00 -6.3518285952025531e-18 -1.3636822777448101e-18 -1.4441625743554549e-02 3.8797851714535434e+00 -3.8802125261354576e-18 -2.7058404544835651e-01 5.0846848186579174e-01 4.0806630199703786e+00 +3.9448633066729739e+00 -3.6344794672422329e-18 1.1581293029516294e-17 +1.6817409381958889e-01 4.0320777239133783e+00 -8.7954588948930866e-19 +1.4551929935812052e-01 3.9989012843738053e-01 3.9977352939318109e+00 C Si 4 4 Cartesian - -0.1284024102 0.0411632429 2.5392886616 - 1.6605207719 -0.0253819250 0.3618554716 - 0.1640583101 2.1563760674 0.4480081454 - 2.1989105171 1.7886495520 2.0950553386 - -0.1078475228 0.2108036732 -0.3228714337 - 1.8447569702 2.1025290971 -0.4494344734 - 1.6511638763 -0.0784401477 1.9721672179 - -0.0382262924 1.6485487817 1.9674381019 + 0.0946327782 -0.2781320996 1.6927633385 + 2.1395921916 -0.5505295607 0.3561009428 + -0.2623111485 2.0403483142 0.5289485573 + 2.0291102621 2.0166734366 1.4189102698 + -0.4274609464 0.4275845699 -0.2870878166 + 1.8398859643 1.9175800522 -0.7344827910 + 1.5794721062 0.2698516856 1.5600233361 + 0.2272461445 2.4170250626 2.4545090907 diff --git a/tests/poscars/POSCAR.SiC.normal b/tests/poscars/POSCAR.SiC.normal index cb047ce5f..e5e2be845 100644 --- a/tests/poscars/POSCAR.SiC.normal +++ b/tests/poscars/POSCAR.SiC.normal @@ -1,16 +1,16 @@ C4 Si4 1.0 -4.0972699470293215e+00 1.0149069314692106e-18 -1.2199095357777910e-17 -1.0494335660165277e-01 4.0313628735769402e+00 -4.8506301604322156e-19 -4.3100124938951284e-01 4.1459063237657107e-01 4.1291342816555661e+00 +4.0354487481064565e+00 1.1027270790560616e-17 2.5642993008475204e-17 +2.0693526054669642e-01 4.1066892997402196e+00 -8.6715682899078028e-18 +4.2891472979598610e-01 5.5796885749827474e-01 4.1100061517204542e+00 C Si 4 4 Cartesian - -0.3125685718 0.3578306188 2.6198865483 - 2.0773182946 -0.1374738727 0.1909202951 - -0.0358416028 2.2809266017 0.5372668691 - 2.1005017067 2.2434342103 2.4385806771 - -0.1565383829 0.1676212247 -0.0494987290 - 1.7296257122 1.9113469908 -0.0379209714 - 2.4121509301 0.3178978635 2.4128037669 - -0.0098316038 1.8287827235 2.0661147519 + 0.2840021179 -0.6003817253 2.7192043932 + 1.7870636884 0.4043686026 -0.3468222936 + -0.4647946470 2.0243215781 -0.0624256251 + 2.1747504767 2.0704389063 2.4228997681 + -0.3817468326 0.0611851697 0.1116535817 + 2.1483270977 2.3528212071 -0.3719335435 + 2.4540854549 0.7297685673 1.8434611305 + 0.1075358778 2.0300985762 2.0687710181 diff --git a/tests/poscars/POSCAR.SiC.uniform b/tests/poscars/POSCAR.SiC.uniform index 7971b3653..9a6214267 100644 --- a/tests/poscars/POSCAR.SiC.uniform +++ b/tests/poscars/POSCAR.SiC.uniform @@ -1,16 +1,16 @@ C4 Si4 1.0 -4.0393044668683302e+00 3.7866840715513190e-18 -3.4832403817178164e-18 -2.4199759881076477e-01 4.1312396844887331e+00 2.2713212548364351e-17 -4.4481025327383367e-01 2.0732581915669496e-01 4.0512745001512140e+00 +4.0817852422279959e+00 -8.0334347485565818e-20 1.4195449766884313e-17 +1.2315711101650186e-01 3.7860044730134366e+00 -1.7522489439993958e-18 +4.0639226354941010e-01 3.2444019947671882e-01 4.2124083393682454e+00 C Si 4 4 Cartesian - -0.2187320173 -0.3278268312 1.9808168543 - 1.7325893560 -0.2907437704 -0.1316407092 - 0.1697204408 2.0456083784 -0.0881243869 - 2.2762667222 2.3246948417 1.8605727611 - -0.0627770396 0.1611339796 0.0659268519 - 2.0190368430 1.8199315498 -0.3333978387 - 2.2452379290 -0.0729729887 2.3597409002 - 0.0376749365 1.7102483986 2.1298706088 + -0.1022885809 0.5553879154 2.3401729917 + 1.4574918226 -0.0391762715 0.0190932943 + -0.0904993945 1.7684440711 -0.1710987880 + 2.0638159391 2.1681985074 2.4580083474 + -0.3261354126 0.1759255681 0.2236316249 + 2.0738703301 2.3024191125 -0.5031567813 + 2.3488610334 0.3825072268 1.7330866249 + -0.0534780340 1.5938442575 1.8486892947 diff --git a/tests/test_perturb.py b/tests/test_perturb.py index 50ca295fa..e34626aab 100644 --- a/tests/test_perturb.py +++ b/tests/test_perturb.py @@ -18,14 +18,14 @@ def rand(self,number): @staticmethod def get_randn_generator(): data = np.asarray([ - [-1.67097394, 1.00565478, 1.75499669], - [ 0.40123822, -0.26832357, 0.51643976], - [-0.70325165, 0.50154408, 1.98677553], - [-0.13296745, 0.55434201, 1.54581177], - [-0.44232103, 0.25854521, -0.43925698], - [-0.86827367, -0.27676856, 0.44365476], - [0.92187855, 1.09601508, 1.56535943], - [-0.61866517, -0.96390108, -0.09593058]]) + [ 0.71878148, -2.20667426, 1.49373955], + [-0.42728113, 1.43836059, -1.17553854], + [-1.70793073, -0.39588759, -0.40880927], + [ 0.17078291, -0.34856352, 1.04307936], + [-0.99103413, -0.1886479, 0.13813131], + [ 0.5839343, 1.04612646, -0.62631026], + [ 0.9752889, 1.85932517, -0.47875828], + [-0.23977172, -0.38373444, -0.04375488]]) count = 0 while True: yield data[count] @@ -33,7 +33,7 @@ def get_randn_generator(): @staticmethod def get_rand_generator(): - yield np.asarray([0.38114487,0.68166018,0.72880982,0.59149848,0.9529526,0.4403747 ]) + yield np.asarray([0.23182233, 0.87106847, 0.68728511, 0.94180274, 0.92860453, 0.69191187]) class UniformGenerator(object): def __init__(self): @@ -46,14 +46,14 @@ def rand(self,number): @staticmethod def get_randn_generator(): - data = [[-0.6338632, -0.80482399, -0.51334222], - [-1.18911393, -1.24521526, -1.34214803], - [ 0.20640311, -0.52666381, -0.71929358], - [ 0.51478689, 1.11296436, -0.61690281], - [-0.77160527, 1.11681373, 0.31440837], - [ 0.71732744, -1.56937797, -1.60701756], - [ 0.6121493, -0.06944649, 2.58409679], - [-0.86858833, -1.44664272, 0.40813879]] + data = [[-0.19313281, 0.80194715, 0.14050915], + [-1.47859926, 0.12921667, -0.17632456], + [-0.60836805, -0.7700423, -0.8386948 ], + [-0.03236753, 0.36690245, 0.5041072 ], + [-1.59366933, 0.37069227, 0.89608291], + [ 0.18165617, 0.53875315, -0.42233955], + [ 0.74052496, 1.26627555, -1.12094823], + [-0.89610092, -1.44247021, -1.3502529 ]] yield np.asarray([0.0001,0.0001,0.0001]) # test for not using small vector count = 0 while True: @@ -62,12 +62,12 @@ def get_randn_generator(): @staticmethod def get_rand_generator(): - data = np.asarray([[0.96216328], [0.50364801], - [0.10566744], [0.82063694], - [0.00212345], [0.3384801], - [0.75329772], [0.59874454]]) + data = np.asarray([[0.71263084], [0.61339295], + [0.22948181], [0.36087632], + [0.17582222], [0.97926742], + [0.84706761], [0.44495513]]) - yield np.asarray([0.24009826,0.9346383,0.54357125,0.05932848,0.96751515,0.78645846]) + yield np.asarray([0.34453551, 0.0618966, 0.9327273, 0.43013654, 0.88624993, 0.48827425]) count =0 while True: yield np.asarray(data[count]) @@ -84,14 +84,14 @@ def rand(self,number): @staticmethod def get_randn_generator(): - data = [[-0.37656207, -0.50545915, 1.23769463], - [-0.77863147, 0.05874498, 0.81611906], - [0.557889, 0.36082455, 1.89695245], - [ 1.07846313, -1.07323247, 0.18513304], - [-0.19418326, 0.75757126, -1.90317556], - [ 0.16438486, 1.40294282, -1.6183854 ], - [-1.71033533, -0.73018772, -0.19946578], - [-0.29332031, -2.00584426, -0.49622181]] + data = np.asarray([[ 0.95410606, -1.62338002, -2.05359934], + [ 0.69213769, -1.26008667, 0.77970721], + [-1.77926476, -0.39227219, 2.31677298], + [ 0.08785233, -0.03966649, -0.45325656], + [-0.53860887, 0.42536802, -0.46167309], + [-0.26865791, -0.19901684, -2.51444768], + [-0.31627314, 0.22076982, -0.36032225], + [0.66731887, 1.2505806, 1.46112938]]) yield np.asarray([0.0001,0.0001,0.0001]) # test for not using small vector count = 0 while True: @@ -100,7 +100,7 @@ def get_randn_generator(): @staticmethod def get_rand_generator(): - yield np.asarray([0.10497343,0.30020419,0.60592304, 0.85344648, 0.60080211, 0.22293347]) + yield np.asarray([0.01525907, 0.68387374, 0.39768541, 0.55596047, 0.26557088, 0.60883073]) # %% class TestPerturbNormal(unittest.TestCase, CompSys): @@ -109,7 +109,7 @@ def setUp (self, random_mock): random_mock.rand = NormalGenerator().rand random_mock.randn = NormalGenerator().randn system_1_origin = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') - self.system_1 = system_1_origin.perturb(1,0.05,0.1,'normal') + self.system_1 = system_1_origin.perturb(1,0.05,0.6,'normal') self.system_2 = dpdata.System('poscars/POSCAR.SiC.normal',fmt='vasp/poscar') self.places = 6 @@ -119,7 +119,7 @@ def setUp (self, random_mock) : random_mock.rand = UniformGenerator().rand random_mock.randn = UniformGenerator().randn system_1_origin = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') - self.system_1 = system_1_origin.perturb(1,0.05,0.1,'uniform') + self.system_1 = system_1_origin.perturb(1,0.05,0.6,'uniform') self.system_2 = dpdata.System('poscars/POSCAR.SiC.uniform',fmt='vasp/poscar') self.places = 6 @@ -129,7 +129,7 @@ def setUp (self, random_mock) : random_mock.rand = ConstGenerator().rand random_mock.randn = ConstGenerator().randn system_1_origin = dpdata.System('poscars/POSCAR.SiC',fmt='vasp/poscar') - self.system_1 = system_1_origin.perturb(1,0.05,0.1,'const') + self.system_1 = system_1_origin.perturb(1,0.05,0.6,'const') self.system_2 = dpdata.System('poscars/POSCAR.SiC.const',fmt='vasp/poscar') self.places = 6 From 185ecf0d01576706682e70272c6a037b0c513e5d Mon Sep 17 00:00:00 2001 From: Han Wang Date: Sat, 14 Dec 2019 17:51:12 +0800 Subject: [PATCH 8/9] Update README.md --- README.md | 29 ++++------------------------- 1 file changed, 4 insertions(+), 25 deletions(-) diff --git a/README.md b/README.md index 85f1f3387..3ec016a47 100644 --- a/README.md +++ b/README.md @@ -153,40 +153,19 @@ dpdata.LabeledSystem('OUTCAR').sub_system([0,-1]).to_deepmd_raw('dpmd_raw') ``` by which only the first and last frames are dumped to `dpmd_raw`. -# replicate +## replicate dpdata will create a super cell of the current atom configuration. - ```python dpdata.System('./POSCAR').replicate((1,2,3,) ) ``` tuple(1,2,3) means don't copy atom configuration in x direction, make 2 copys in y direction, make 3 copys in z direction. -# perturb -dpdata will disturb the box size and shape and change atom coordinates randomly. +## perturb +By the following example, each frame of the original system (`dpdata.System('./POSCAR')`) is perturbed to generate three new frames. For each frame, the cell is perturbed by 5% and the atom positions are perturbed by 0.6 Angstrom. `atom_pert_style` indicates that the perturbation to the atom positions is subject to normal distribution. Other available options to `atom_pert_style` are`uniform` (uniform in a ball), and `const` (uniform on a sphere). ```python perturbed_system = dpdata.System('./POSCAR').perturb(pert_num=3, - box_pert_fraction=0.05, + cell_pert_fraction=0.05, atom_pert_distance=0.6, atom_pert_style='normal') print(perturbed_system.data) ``` - -### `pert_num` - -Each frame in the input system will generate `pert_num` frames. -That means the command will return a system containing `frames of input system * pert_num` frames. - -### `box_pert_fraction` -A relative length that determines the length the box size change in each frame. It is just a fraction and doesn't have unit. Typlicaly, for cubic box with side length `side_length` , `side_length` will increase or decrease a random value from the intervel `[-side_length*box_pert_fraction, side_length*box_pert_fraction]`. - -It will also change the shape of the box. That means an orthogonal box will become a non-orthogonal box after perturbing. The angle of inclination of the box is a random variable. -`box_pert_fraction` is also relating to the probability distribution function of the angle. - -### `atom_pert_distance` -unit:Angstrom. Determine the distance atoms move in each frame. -The mean value of the distance that atom moves is about `atom_pert_distance`. - -### `atom_pert_style` -The probability distribution function used to change atom coordinates. -available options:`'uniform', 'normal', 'const'`. -The direction atoms move and box deformation is random. From 0e4d5c342c259de5c1d3e6c841156c0739c1b3a4 Mon Sep 17 00:00:00 2001 From: Han Wang Date: Sat, 14 Dec 2019 18:11:50 +0800 Subject: [PATCH 9/9] change all 'box' to 'cell' --- dpdata/system.py | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/dpdata/system.py b/dpdata/system.py index ab020abc5..ebc38ea1b 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -626,13 +626,12 @@ def replicate(self, ncopy): def perturb(self, pert_num, - box_pert_fraction, + cell_pert_fraction, atom_pert_distance, atom_pert_style='normal'): """ Perturb each frame in the system randomly. - The box will be changed randomly, - and atoms will move random distance in random direction. + The cell will be deformed randomly, and atoms will be displaced by a random distance in random direction. Parameters ---------- @@ -640,12 +639,11 @@ def perturb(self, Each frame in the system will make `pert_num` copies, and all the copies will be perturbed. That means the system to be returned will contain `pert_num` * frame_num of the input system. - box_pert_fraction : float - A fraction determines how much will box deform, typically less than 0.3. - For a cubic box with side length `side_length`, - `side_length` will increase or decrease with the max value `box_pert_fraction*side_length` . - If box_pert_fraction is not zero, the shape of the box will also be changed, - and that means a orthogonal box will become a non-orthogonal one. + cell_pert_fraction : float + A fraction determines how much (relatively) will cell deform. + The cell of each frame is deformed by a symmetric matrix perturbed from identity. + The perturbation to the diagonal part is subject to a uniform distribution in [-cell_pert_fraction, cell_pert_fraction), + and the perturbation to the off-diagonal part is subject to a uniform distribution in [-0.5*cell_pert_fraction, 0.5*cell_pert_fraction). atom_pert_distance: float unit: Angstrom. A distance determines how far atoms will move. Atoms will move about `atom_pert_distance` in random direction. @@ -670,9 +668,9 @@ def perturb(self, for ii in range(nframes): for jj in range(pert_num): tmp_system = self[ii].copy() - box_perturb_matrix = get_box_perturb_matrix(box_pert_fraction) - tmp_system.data['cells'][0] = np.matmul(tmp_system.data['cells'][0],box_perturb_matrix) - tmp_system.data['coords'][0] = np.matmul(tmp_system.data['coords'][0],box_perturb_matrix) + cell_perturb_matrix = get_cell_perturb_matrix(cell_pert_fraction) + tmp_system.data['cells'][0] = np.matmul(tmp_system.data['cells'][0],cell_perturb_matrix) + tmp_system.data['coords'][0] = np.matmul(tmp_system.data['coords'][0],cell_perturb_matrix) for kk in range(len(tmp_system.data['coords'][0])): atom_perturb_vector = get_atom_perturb_vector(atom_pert_distance, atom_pert_style) tmp_system.data['coords'][0][kk] += atom_perturb_vector @@ -680,17 +678,17 @@ def perturb(self, perturbed_system.append(tmp_system) return perturbed_system -def get_box_perturb_matrix(box_pert_fraction): - if box_pert_fraction<0: - raise RuntimeError('box_pert_fraction can not be negative') +def get_cell_perturb_matrix(cell_pert_fraction): + if cell_pert_fraction<0: + raise RuntimeError('cell_pert_fraction can not be negative') e0 = np.random.rand(6) - e = e0 * 2 *box_pert_fraction - box_pert_fraction - box_pert_matrix = np.array( + e = e0 * 2 *cell_pert_fraction - cell_pert_fraction + cell_pert_matrix = np.array( [[1+e[0], 0.5 * e[5], 0.5 * e[4]], [0.5 * e[5], 1+e[1], 0.5 * e[3]], [0.5 * e[4], 0.5 * e[3], 1+e[2]]] ) - return box_pert_matrix + return cell_pert_matrix def get_atom_perturb_vector(atom_pert_distance, atom_pert_style='normal'): random_vector = None