scnc
/
phonopy
mirror of https://github.com/phonopy/phonopy.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

C implementation of permulation symmetry fc4

This commit is contained in:
Atsushi Togo 2013-08-28 16:22:36 +09:00
parent c8f21568ff
commit 3608972dc3
7 changed files with 393 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
anharmonic/phonon3/fc3.py
View File

@ -1,6 +1,6 @@
import sys
import numpy as np
from phonopy.harmonic.force_constants import similarity_transformation, set_permutation_symmetry, set_translational_invariance_per_index, distribute_force_constants, solve_force_constants, get_rotated_displacement, get_positions_sent_by_rot_inv
from phonopy.harmonic.force_constants import similarity_transformation, set_permutation_symmetry, distribute_force_constants, solve_force_constants, get_rotated_displacement, get_positions_sent_by_rot_inv, set_translational_invariance
from anharmonic.phonon3.displacement_fc3 import get_reduced_site_symmetry, get_bond_symmetry
from anharmonic.file_IO import write_fc2_dat
@ -287,7 +287,7 @@ def get_constrained_fc2(supercell,
symprec)
if is_translational_symmetry:
set_translational_invariance_per_index(fc2)
set_translational_invariance(fc2)
if is_permutation_symmetry:
set_permutation_symmetry(fc2)

90
anharmonic/phonon4/fc4.py
View File

@ -89,15 +89,20 @@ def set_translational_invariance_fc4_per_index_py(fc4, index=0):
fc4[i, j, k, :, l, m, n, p]) / fc4.shape[3]
def set_permutation_symmetry_fc4(fc4):
num_atom = fc4.shape[0]
fc4_sym = np.zeros(fc4.shape, dtype='double')
for (i, j, k, l) in list(
np.ndindex(num_atom, num_atom, num_atom, num_atom)):
fc4_sym[i, j, k, l] = set_permutation_symmetry_fc4_part(fc4, i, j, k, l)
for (i, j, k, l) in list(
np.ndindex(num_atom, num_atom, num_atom, num_atom)):
fc4[i, j, k, l] = fc4_sym[i, j, k, l]
try:
import anharmonic._phono4py as phono4c
phono4c.permutation_symmetry_fc4(fc4)
except ImportError:
num_atom = fc4.shape[0]
fc4_sym = np.zeros(fc4.shape, dtype='double')
for (i, j, k, l) in list(
np.ndindex(num_atom, num_atom, num_atom, num_atom)):
fc4_sym[i, j, k, l] = set_permutation_symmetry_fc4_part(
fc4, i, j, k, l)
for (i, j, k, l) in list(
np.ndindex(num_atom, num_atom, num_atom, num_atom)):
fc4[i, j, k, l] = fc4_sym[i, j, k, l]
def set_permutation_symmetry_fc4_part(fc4, a, b, c, d):
tensor4 = np.zeros((3, 3, 3, 3), dtype='double')
@ -131,25 +136,33 @@ def set_permutation_symmetry_fc4_part(fc4, a, b, c, d):
return tensor4
def show_drift_fc4(fc4, name="fc4"):
num_atom = fc4.shape[0]
maxval1 = 0
maxval2 = 0
maxval3 = 0
maxval4 = 0
for i, j, k, l, m, n, p in list(
np.ndindex((num_atom, num_atom, num_atom, 3, 3, 3, 3))):
val1 = fc4[:, i, j, k, l, m, n, p].sum()
val2 = fc4[k, :, i, j, l, m, n, p].sum()
val3 = fc4[j, k, :, i, l, m, n, p].sum()
val4 = fc4[i, j, k, :, l, m, n, p].sum()
if abs(val1) > abs(maxval1):
maxval1 = val1
if abs(val2) > abs(maxval2):
maxval2 = val2
if abs(val3) > abs(maxval3):
maxval3 = val3
if abs(val4) > abs(maxval4):
maxval4 = val4
try:
import anharmonic._phono4py as phono4c
(maxval1,
maxval2,
maxval3,
maxval4) = phono4c.drift_fc4(fc4)
except ImportError:
num_atom = fc4.shape[0]
maxval1 = 0
maxval2 = 0
maxval3 = 0
maxval4 = 0
for i, j, k, l, m, n, p in list(
np.ndindex((num_atom, num_atom, num_atom, 3, 3, 3, 3))):
val1 = fc4[:, i, j, k, l, m, n, p].sum()
val2 = fc4[k, :, i, j, l, m, n, p].sum()
val3 = fc4[j, k, :, i, l, m, n, p].sum()
val4 = fc4[i, j, k, :, l, m, n, p].sum()
if abs(val1) > abs(maxval1):
maxval1 = val1
if abs(val2) > abs(maxval2):
maxval2 = val2
if abs(val3) > abs(maxval3):
maxval3 = val3
if abs(val4) > abs(maxval4):
maxval4 = val4
print ("max drift of %s:" % name), maxval1, maxval2, maxval3, maxval4
def distribute_fc4(fc4,
@ -310,7 +323,8 @@ def _get_delta_fc3(dataset_first_atom,
symprec,
verbose)
show_drift_fc3(disp_fc3, name="fc3 with disp.")
if verbose:
show_drift_fc3(disp_fc3, name="fc3 with disp.")
return disp_fc3 - fc3
@ -396,6 +410,7 @@ def _get_constrained_fc3(supercell,
(atom1 + 1, atom2 + 1))
for i, v in enumerate(disps2):
print " [%7.4f %7.4f %7.4f]" % tuple(v)
solve_fc3(fc3,
atom2,
supercell,
@ -449,24 +464,25 @@ def _solve_fc4(fc4,
rot_disps = get_rotated_displacement(displacements_first, site_sym_cart)
inv_U = np.linalg.pinv(rot_disps)
for (i, j, k) in list(np.ndindex(num_atom, num_atom, num_atom)):
fc4[first_atom_num, i, j, k] = np.dot(
inv_U, _rotate_delta_fc3s(
i, j, k, delta_fc3s, rot_map_syms, site_sym_cart)
).reshape(3, 3, 3, 3)
def _rotate_delta_fc3s(i, j, k, delta_fc3s, rot_map_syms, site_sym_cart):
rotated_fc3s = np.zeros((len(delta_fc3s), len(site_sym_cart), 3, 3, 3),
dtype='double')
try:
import anharmonic._phono4py as phono4c
phono4c.rotate_delta_fc3s(rotated_fc3s,
delta_fc3s,
rot_map_syms,
site_sym_cart,
i,
j,
k)
phono4c.rotate_delta_fc3s_elem(rotated_fc3s,
delta_fc3s,
rot_map_syms,
site_sym_cart,
i,
j,
k)
return np.reshape(rotated_fc3s, (-1, 27))
except ImportError:
print "Copying delta fc3s at (%d, %d, %d)" % (i + 1, j + 1, k + 1)

74
c/_phono4py.c
View File

@ -14,9 +14,12 @@ static PyObject * py_real_to_reciprocal4(PyObject *self, PyObject *args);
static PyObject * py_reciprocal_to_normal4(PyObject *self, PyObject *args);
static PyObject * py_set_phonons_grid_points(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc4(PyObject *self, PyObject *args);
static PyObject * py_rotate_delta_fc3s(PyObject *self, PyObject *args);
static PyObject * py_rotate_delta_fc3s_elem(PyObject *self, PyObject *args);
static PyObject * py_set_translational_invariance_fc4(PyObject *self,
PyObject *args);
static PyObject * py_set_permutation_symmetry_fc4(PyObject *self,
PyObject *args);
static PyObject * py_get_drift_fc4(PyObject *self, PyObject *args);
static PyMethodDef functions[] = {
{"fc4_normal_for_frequency_shift", py_get_fc4_normal_for_frequency_shift, METH_VARARGS, "Calculate fc4 normal for frequency shift"},
@ -25,8 +28,10 @@ static PyMethodDef functions[] = {
{"reciprocal_to_normal4", py_reciprocal_to_normal4, METH_VARARGS, "Transform fc4 of reciprocal space to normal coordinate in special case for frequency shift"},
{"phonons_grid_points", py_set_phonons_grid_points, METH_VARARGS, "Set phonons on grid points"},
{"distribute_fc4", py_distribute_fc4, METH_VARARGS, "Distribute least fc4 to full fc4"},
{"rotate_delta_fc3s", py_rotate_delta_fc3s, METH_VARARGS, "Rotate delta fc3s"},
{"rotate_delta_fc3s_elem", py_rotate_delta_fc3s_elem, METH_VARARGS, "Rotate delta fc3s for a set of atomic indices"},
{"translational_invariance_fc4", py_set_translational_invariance_fc4, METH_VARARGS, "Set translational invariance for fc4"},
{"permutation_symmetry_fc4", py_set_permutation_symmetry_fc4, METH_VARARGS, "Set permutation symmetry for fc4"},
{"drift_fc4", py_get_drift_fc4, METH_VARARGS, "Get drifts of fc4"},
{NULL, NULL, 0, NULL}
};
@ -392,7 +397,7 @@ static PyObject * py_distribute_fc4(PyObject *self, PyObject *args)
rot_cart_inv));
}
static PyObject * py_rotate_delta_fc3s(PyObject *self, PyObject *args)
static PyObject * py_rotate_delta_fc3s_elem(PyObject *self, PyObject *args)
{
PyArrayObject* rotated_delta_fc3s_py;
PyArrayObject* delta_fc3s_py;
@ -419,16 +424,16 @@ static PyObject * py_rotate_delta_fc3s(PyObject *self, PyObject *args)
const int num_delta_fc3s = (int)delta_fc3s_py->dimensions[0];
const int num_atom = (int)delta_fc3s_py->dimensions[1];
return PyInt_FromLong((long) rotate_delta_fc3s(rotated_delta_fc3s,
delta_fc3s,
rot_map_syms,
site_syms_cart,
num_rot,
num_delta_fc3s,
atom1,
atom2,
atom3,
num_atom));
return PyInt_FromLong((long) rotate_delta_fc3s_elem(rotated_delta_fc3s,
delta_fc3s,
rot_map_syms,
site_syms_cart,
num_rot,
num_delta_fc3s,
atom1,
atom2,
atom3,
num_atom));
}
static PyObject * py_set_translational_invariance_fc4(PyObject *self,
@ -450,3 +455,46 @@ static PyObject * py_set_translational_invariance_fc4(PyObject *self,
Py_RETURN_NONE;
}
static PyObject * py_set_permutation_symmetry_fc4(PyObject *self, PyObject *args)
{
PyArrayObject* fc4_py;
if (!PyArg_ParseTuple(args, "O",
&fc4_py)) {
return NULL;
}
double* fc4 = (double*)fc4_py->data;
const int num_atom = (int)fc4_py->dimensions[0];
set_permutation_symmetry_fc4(fc4, num_atom);
Py_RETURN_NONE;
}
static PyObject * py_get_drift_fc4(PyObject *self, PyObject *args)
{
PyArrayObject* fc4_py;
if (!PyArg_ParseTuple(args, "O",
&fc4_py)) {
return NULL;
}
double* fc4 = (double*)fc4_py->data;
const int num_atom = (int)fc4_py->dimensions[0];
int i;
double drift[4];
PyObject* drift_py;
get_drift_fc4(drift, fc4, num_atom);
drift_py = PyList_New(4);
for (i = 0; i < 4; i++) {
PyList_SetItem(drift_py, i, PyFloat_FromDouble(drift[i]));
}
return drift_py;
}

1
c/anharmonic/phonon3/fc3.c
View File

@ -17,6 +17,7 @@ int distribute_fc3(double *fc3,
third_atom_rot = atom_mapping[third_atom];
#pragma omp parallel for private(j, atom_rot_i, atom_rot_j, tensor)
for (i = 0; i < num_atom; i++) {
atom_rot_i = atom_mapping[i];

305
c/anharmonic/phonon4/fc4.c
View File

@ -1,3 +1,5 @@
#include <math.h>
#include <stdlib.h>
#include "phonon3_h/fc3.h"
#include "phonon4_h/fc4.h"
@ -19,17 +21,31 @@ static double tensor4_rotation_elem(const double tensor[81],
const int n,
const int p,
const int q);
static double get_drift_fc4_elem(const double *fc4,
const int num_atom,
const int i,
const int j,
const int k,
const int l,
const int index);
void set_permutation_symmetry_fc4_elem(double *fc4_elem,
const double *fc4_tmp,
const int a,
const int b,
const int c,
const int d,
const int num_atom);
int rotate_delta_fc3s(double *rotated_delta_fc3s,
const double *delta_fc3s,
const int *rot_map_syms,
const double *site_sym_cart,
const int num_rot,
const int num_delta_fc3s,
const int atom1,
const int atom2,
const int atom3,
const int num_atom)
int rotate_delta_fc3s_elem(double *rotated_delta_fc3s,
const double *delta_fc3s,
const int *rot_map_syms,
const double *site_sym_cart,
const int num_rot,
const int num_delta_fc3s,
const int atom1,
const int atom2,
const int atom3,
const int num_atom)
{
int i, j;
double *rot_tensor;
@ -64,7 +80,7 @@ int distribute_fc4(double *fc4,
fourth_atom_rot = atom_mapping[fourth_atom];
#pragma omp parallel for private(j, k, atom_rot_i, atom_rot_j, atom_rot_k)
#pragma omp parallel for private(j, k, atom_rot_i, atom_rot_j, atom_rot_k, tensor)
for (i = 0; i < num_atom; i++) {
atom_rot_i = atom_mapping[i];
@ -112,44 +128,15 @@ void set_translational_invariance_fc4_per_index(double *fc4,
const int index)
{
int i, j, k, l, m;
double sum, drift;
double drift;
#pragma omp parallel for private(j, k, l, m, sum, drift)
#pragma omp parallel for private(j, k, l, m, drift)
for (i = 0; i < 81; i++) {
for (j = 0; j < num_atom; j++) {
for (k = 0; k < num_atom; k++) {
for (l = 0; l < num_atom; l++) {
sum = 0;
for (m = 0; m < num_atom; m++) {
switch (index) {
case 0:
sum += fc4[m * num_atom * num_atom * num_atom * 81 +
j * num_atom * num_atom * 81 +
k * num_atom * 81 +
l * 81 + i];
break;
case 1:
sum += fc4[j * num_atom * num_atom * num_atom * 81 +
m * num_atom * num_atom * 81 +
k * num_atom * 81 +
l * 81 + i];
break;
case 2:
sum += fc4[j * num_atom * num_atom * num_atom * 81 +
k * num_atom * num_atom * 81 +
m * num_atom * 81 +
l * 81 + i];
break;
case 3:
sum += fc4[j * num_atom * num_atom * num_atom * 81 +
k * num_atom * num_atom * 81 +
l * num_atom * 81 +
m * 81 + i];
break;
}
}
drift = sum / num_atom;
drift =
get_drift_fc4_elem(fc4, num_atom, i, j, k, l, index) / num_atom;
for (m = 0; m < num_atom; m++) {
switch (index) {
case 0:
@ -182,9 +169,237 @@ void set_translational_invariance_fc4_per_index(double *fc4,
}
}
}
}
void set_permutation_symmetry_fc4(double *fc4, const int num_atom)
{
double *fc4_tmp;
int i, j, k, l;
fc4_tmp = (double*)malloc(sizeof(double) *
num_atom * num_atom * num_atom * num_atom * 81);
#pragma omp parallel for
for (i = 0; i < num_atom * num_atom * num_atom * num_atom * 81; i++) {
fc4_tmp[i] = fc4[i];
}
#pragma omp parallel for private(j, k, l)
for (i = 0; i < num_atom; i++) {
for (j = 0; j < num_atom; j++) {
for (k = 0; k < num_atom; k++) {
for (l = 0; l < num_atom; l++) {
set_permutation_symmetry_fc4_elem
(fc4 +
i * num_atom * num_atom * num_atom * 81 +
j * num_atom * num_atom * 81 +
k * num_atom * 81 +
l * 81, fc4_tmp, i, j, k, l, num_atom);
}
}
}
}
free(fc4_tmp);
}
void set_permutation_symmetry_fc4_elem(double *fc4_elem,
const double *fc4_tmp,
const int a,
const int b,
const int c,
const int d,
const int num_atom)
{
int i, j, k, l;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
fc4_elem[i * 27 + j * 9 + k * 3 + l] =
(fc4_tmp[a * num_atom * num_atom * num_atom * 81 +
b * num_atom * num_atom * 81 +
c * num_atom * 81+
d * 81 + i * 27 + j * 9 + k * 3 + l] +
fc4_tmp[a * num_atom * num_atom * num_atom * 81 +
b * num_atom * num_atom * 81 +
d * num_atom * 81+
c * 81 + i * 27 + j * 9 + l * 3 + k] +
fc4_tmp[a * num_atom * num_atom * num_atom * 81 +
c * num_atom * num_atom * 81 +
b * num_atom * 81+
d * 81 + i * 27 + k * 9 + j * 3 + l] +
fc4_tmp[a * num_atom * num_atom * num_atom * 81 +
c * num_atom * num_atom * 81 +
d * num_atom * 81+
b * 81 + i * 27 + k * 9 + l * 3 + j] +
fc4_tmp[a * num_atom * num_atom * num_atom * 81 +
d * num_atom * num_atom * 81 +
b * num_atom * 81+
c * 81 + i * 27 + l * 9 + j * 3 + k] +
fc4_tmp[a * num_atom * num_atom * num_atom * 81 +
d * num_atom * num_atom * 81 +
c * num_atom * 81+
b * 81 + i * 27 + l * 9 + k * 3 + j] +
fc4_tmp[b * num_atom * num_atom * num_atom * 81 +
a * num_atom * num_atom * 81 +
c * num_atom * 81+
d * 81 + j * 27 + i * 9 + k * 3 + l] +
fc4_tmp[b * num_atom * num_atom * num_atom * 81 +
a * num_atom * num_atom * 81 +
d * num_atom * 81+
c * 81 + j * 27 + i * 9 + l * 3 + k] +
fc4_tmp[b * num_atom * num_atom * num_atom * 81 +
c * num_atom * num_atom * 81 +
a * num_atom * 81+
d * 81 + j * 27 + k * 9 + i * 3 + l] +
fc4_tmp[b * num_atom * num_atom * num_atom * 81 +
c * num_atom * num_atom * 81 +
d * num_atom * 81+
a * 81 + j * 27 + k * 9 + l * 3 + i] +
fc4_tmp[b * num_atom * num_atom * num_atom * 81 +
d * num_atom * num_atom * 81 +
a * num_atom * 81+
c * 81 + j * 27 + l * 9 + i * 3 + k] +
fc4_tmp[b * num_atom * num_atom * num_atom * 81 +
d * num_atom * num_atom * 81 +
c * num_atom * 81+
a * 81 + j * 27 + l * 9 + k * 3 + i] +
fc4_tmp[c * num_atom * num_atom * num_atom * 81 +
a * num_atom * num_atom * 81 +
b * num_atom * 81+
d * 81 + k * 27 + i * 9 + j * 3 + l] +
fc4_tmp[c * num_atom * num_atom * num_atom * 81 +
a * num_atom * num_atom * 81 +
d * num_atom * 81+
b * 81 + k * 27 + i * 9 + l * 3 + j] +
fc4_tmp[c * num_atom * num_atom * num_atom * 81 +
b * num_atom * num_atom * 81 +
a * num_atom * 81+
d * 81 + k * 27 + j * 9 + i * 3 + l] +
fc4_tmp[c * num_atom * num_atom * num_atom * 81 +
b * num_atom * num_atom * 81 +
d * num_atom * 81+
a * 81 + k * 27 + j * 9 + l * 3 + i] +
fc4_tmp[c * num_atom * num_atom * num_atom * 81 +
d * num_atom * num_atom * 81 +
a * num_atom * 81+
b * 81 + k * 27 + l * 9 + i * 3 + j] +
fc4_tmp[c * num_atom * num_atom * num_atom * 81 +
d * num_atom * num_atom * 81 +
b * num_atom * 81+
a * 81 + k * 27 + l * 9 + j * 3 + i] +
fc4_tmp[d * num_atom * num_atom * num_atom * 81 +
a * num_atom * num_atom * 81 +
b * num_atom * 81+
c * 81 + l * 27 + i * 9 + j * 3 + k] +
fc4_tmp[d * num_atom * num_atom * num_atom * 81 +
a * num_atom * num_atom * 81 +
c * num_atom * 81+
b * 81 + l * 27 + i * 9 + k * 3 + j] +
fc4_tmp[d * num_atom * num_atom * num_atom * 81 +
b * num_atom * num_atom * 81 +
a * num_atom * 81+
c * 81 + l * 27 + j * 9 + i * 3 + k] +
fc4_tmp[d * num_atom * num_atom * num_atom * 81 +
b * num_atom * num_atom * 81 +
c * num_atom * 81+
a * 81 + l * 27 + j * 9 + k * 3 + i] +
fc4_tmp[d * num_atom * num_atom * num_atom * 81 +
c * num_atom * num_atom * 81 +
a * num_atom * 81+
b * 81 + l * 27 + k * 9 + i * 3 + j] +
fc4_tmp[d * num_atom * num_atom * num_atom * 81 +
c * num_atom * num_atom * 81 +
b * num_atom * 81+
a * 81 + l * 27 + k * 9 + j * 3 + i]) / 24;
}
}
}
}
}
void get_drift_fc4(double *drifts_out, const double *fc4, const int num_atom)
{
int i, j, k, l, index;
double drift;
double max_drift[81];
for (index = 0; index < 4; index++) {
for (i = 0; i < 81; i++) {
max_drift[i] = 0;
}
#pragma omp parallel for private(j, k, l, drift)
for (i = 0; i < 81; i++) {
for (j = 0; j < num_atom; j++) {
for (k = 0; k < num_atom; k++) {
for (l = 0; l < num_atom; l++) {
drift = get_drift_fc4_elem(fc4, num_atom, i, j, k, l, index);
if (fabs(max_drift[i]) < fabs(drift)) {
max_drift[i] = drift;
}
}
}
}
}
drift = 0;
for (i = 0; i < 81; i++) {
if (fabs(drift) < fabs(max_drift[i])) {
drift = max_drift[i];
}
}
drifts_out[index] = drift;
}
}
static double get_drift_fc4_elem(const double *fc4,
const int num_atom,
const int i,
const int j,
const int k,
const int l,
const int index)
{
double sum;
int m;
sum = 0;
for (m = 0; m < num_atom; m++) {
switch (index) {
case 0:
sum += fc4[m * num_atom * num_atom * num_atom * 81 +
j * num_atom * num_atom * 81 +
k * num_atom * 81 +
l * 81 + i];
break;
case 1:
sum += fc4[j * num_atom * num_atom * num_atom * 81 +
m * num_atom * num_atom * 81 +
k * num_atom * 81 +
l * 81 + i];
break;
case 2:
sum += fc4[j * num_atom * num_atom * num_atom * 81 +
k * num_atom * num_atom * 81 +
m * num_atom * 81 +
l * 81 + i];
break;
case 3:
sum += fc4[j * num_atom * num_atom * num_atom * 81 +
k * num_atom * num_atom * 81 +
l * num_atom * 81 +
m * 81 + i];
break;
}
}
return sum;
}
static void tensor4_roation(double *rot_tensor,
const double *fc4,
const int atom_i,

23
c/anharmonic_h/phonon4_h/fc4.h
View File

@ -1,16 +1,16 @@
#ifndef __fc4_H__
#define __fc4_H__
int rotate_delta_fc3s(double *rotated_delta_fc3s,
const double *delta_fc3s,
const int *rot_map_syms,
const double *site_sym_cart,
const int num_rot,
const int num_delta_fc3s,
const int atom1,
const int atom2,
const int atom3,
const int num_atom);
int rotate_delta_fc3s_elem(double *rotated_delta_fc3s,
const double *delta_fc3s,
const int *rot_map_syms,
const double *site_sym_cart,
const int num_rot,
const int num_delta_fc3s,
const int atom1,
const int atom2,
const int atom3,
const int num_atom);
int distribute_fc4(double *fc4,
const int fourth_atom,
const int *atom_mapping,
@ -21,5 +21,8 @@ void set_translational_invariance_fc4(double *fc4,
void set_translational_invariance_fc4_per_index(double *fc4,
const int num_atom,
const int index);
void set_permutation_symmetry_fc4(double *fc4, const int num_atom);
void get_drift_fc4(double *drift, const double *fc4, const int num_atom);
#endif

3
scripts/phono4py
View File

@ -403,6 +403,9 @@ else:
print "----- Symmetrize fc4 by index exchange in real space -----"
set_permutation_symmetry_fc4(fc4)
if log_level:
print "(Calculating fc4 drift...)"
show_drift_fc4(fc4)
if not options.read_fc4: