commit 86112fec2aa5fdbdc1ff189b3e32b84878d41459
Author: Marti Aliod Carles <cmartial@pxe.cbp.ens-lyon.fr>
Date:   Wed Oct 2 16:50:18 2019 +0200

    Initial files

diff --git a/DockOnSurf_path b/DockOnSurf_path
new file mode 100644
index 0000000..f0b5880
--- /dev/null
+++ b/DockOnSurf_path
@@ -0,0 +1,33 @@
+#####################################################
+# Please change here the paths to your directories: #
+#####################################################
+
+
+# Directory of the DockOnSurf script:
+export DockOnSurf_path="/home/sblanck/DockOnSurf_hematite"
+
+# Directory where you want to get the results for the molecule alone:
+# The script will automatically create a subdirectory with the name of the studied molecule
+export Molecule_results_path="/home/sblanck/cp2k/hematite_molecules"
+
+# Directory where you want to get the results for the asorption structures:
+# The script will automatically create a subdirectory with the name of the studied molecule
+export MolOnSurf_results_path="/home/sblanck/cp2k/hematite"
+
+# Directory where your surface xyz files are stored:
+export Surface_path="/home/sblanck/DockOnSurf_hematite/surfaces"
+
+# Directory where your molecule mol files are stored:
+export Molecule_path="/home/sblanck/DockOnSurf_hematite/molecules"
+
+# Input file for CP2K calculation for molecule alone:
+# Please check that the PROJECT_NAME is "molecule"
+export CP2K_input_molecule="/home/sblanck/DockOnSurf_hematite/CP2K_input/molecule.inp"
+
+# Input file for CP2K calculation for molecule on the surface:
+# Please check that the PROJECT_NAME is "surf_molecule"
+export CP2K_input_MolOnSurf="/home/sblanck/DockOnSurf_hematite/CP2K_input/adsorption.inp"
+
+# Submission script file for CP2K calculations:
+# Please check that the name of the Job is "XXXX"
+export CP2K_sub="/home/sblanck/DockOnSurf_hematite/CP2K_input/gammakpoint.j"
diff --git a/global_script.sh b/global_script.sh
new file mode 100755
index 0000000..c19e2ac
--- /dev/null
+++ b/global_script.sh
@@ -0,0 +1,178 @@
+#!/bin/bash
+
+user="$(pwd | awk -F"/" '{print $3}')"
+
+file_path="$(find /home/${user} -name "DockOnSurf_path")"
+
+### Mise en place de l'arborescence
+
+source ${file_path}
+
+### Récupération des vecteurs de cellule
+
+a1="$(grep "A\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^A/ {print}' | awk '{print $2}')"
+if [[ $a1 =~ "E+00" ]]; then a1="$(echo $a1 | cut -d "E" -f1)"; fi ; export a1
+a2="$(grep "A\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^A/ {print}' | awk '{print $3}')"
+if [[ $a2 =~ "E+00" ]]; then a2="$(echo $a2 | cut -d "E" -f1)"; fi ; export a2
+a3="$(grep "A\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^A/ {print}' | awk '{print $4}')"
+if [[ $a3 =~ "E+00" ]]; then a3="$(echo $a3 | cut -d "E" -f1)"; fi ; export a3
+b1="$(grep "B\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^B/ {print}' | awk '{print $2}')"
+if [[ $b1 =~ "E+00" ]]; then b1="$(echo $b1 | cut -d "E" -f1)"; fi ; export b1
+b2="$(grep "B\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^B/ {print}' | awk '{print $3}')"
+if [[ $b2 =~ "E+00" ]]; then b2="$(echo $b2 | cut -d "E" -f1)"; fi ; export b2
+b3="$(grep "B\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^B/ {print}' | awk '{print $4}')"
+if [[ $b3 =~ "E+00" ]]; then b3="$(echo $b3 | cut -d "E" -f1)"; fi ; export b3
+c1="$(grep "C\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^C/ {print}' | awk '{print $2}')"
+if [[ $c1 =~ "E+00" ]]; then c1="$(echo $c1 | cut -d "E" -f1)"; fi ; export c1
+c2="$(grep "C\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^C/ {print}' | awk '{print $3}')"
+if [[ $c2 =~ "E+00" ]]; then c2="$(echo $c2 | cut -d "E" -f1)"; fi ; export c2
+c3="$(grep "C\ " ${CP2K_input_MolOnSurf} | awk '$1 ~/^C/ {print}' | awk '{print $4}')"
+if [[ $c3 =~ "E+00" ]]; then c3="$(echo $c3 | cut -d "E" -f1)"; fi ; export c3
+
+### Paramètres initiaux
+
+echo 'Name of the molecule (without .mol) '
+read molecule
+
+echo '' 
+
+module load rdkit
+
+### Generation des conformeres
+
+${DockOnSurf_path}/modules/generation_conformeres.py $molecule 5000
+echo ' -- CONFORMERS HAVE BEEN GENERATED --'
+
+### Conversion des fichiers .mol en .xyz
+
+${DockOnSurf_path}/modules/mol_to_xyz.sh $molecule
+
+### Lancement des calculs pour les molecules seules
+
+${DockOnSurf_path}/modules/launch_cp2k_molecule_seule.sh $molecule
+
+### Attente jusqu'à ce que tous les calculs aient fini
+
+mol="$(echo $molecule | cut -c1-10)"
+while [ "$(qstat | grep $mol | wc -l)" != 0 ];
+do 
+	sleep 300
+done
+
+echo ' -- ALL CALCULATIONS FOR MOLECULE ALONE HAVE BEEN DONE --'
+
+### Suppression des fichiers inutiles
+
+rm -r ${Molecule_results_path}/${molecule}/*/*RESTART*
+rm -r ${Molecule_results_path}/${molecule}/*/*restart*
+rm -r ${Molecule_results_path}/${molecule}/*/*BFGS*
+rm -r ${Molecule_results_path}/${molecule}/*/*.p*
+
+echo ' '
+echo ' -- RESULTS --'
+
+min_energy_molecule=0
+for struct in ${Molecule_results_path}/${molecule}/${molecule}*; do
+	energy_molecule="$(awk '/ENERGY/{E=$NF} END{print E}' ${struct}/*${molecule}.out)"
+	if [[ "${energy_molecule}" > "${min_energy_molecule}" ]]
+	then 
+		min_energy_molecule=${energy_molecule}
+		numero_mol_stable="$(echo $struct | awk -F'/' '{print $NF}')"
+	fi
+done
+echo "   Most stable structure for molecule alone: ${numero_mol_stable} with energy ${min_energy_molecule} Ha"
+
+echo ' ' 
+echo ' -- ARGUMENTS FOR ADSORPTION -- '
+echo ' '
+echo 'Name of the file to use for the surface (without .xyz) '
+read surface
+echo 'Number of the central atom of the molecule'
+read center
+read -a list_atom_1_mol -p 'List of atoms of the molecule to adsorb '
+read -a list_atom_2_mol -p 'List of neighboring atoms of the molecule (linked to the atoms to adsorb) '
+nombre_atom_mol=${#list_atom_1_mol[@]}
+if [ ${#list_atom_2_mol[@]} != ${nombre_atom_mol} ]
+then
+	echo "  Problem: number of atoms different in the two lists"
+	exit
+fi
+read -a list_atom_surf -p 'List of the atoms of the surface that can be adsorption sites '
+echo 'Value of cutoff energy for reoptimisation (eV)'
+read cutoff
+echo '' 
+
+### Adsorption de la molecule sur la surface
+
+${DockOnSurf_path}/modules/script_add_adsorbate+diss.sh $molecule $surface "${list_atom_1_mol[*]}" "${list_atom_2_mol[*]}" "${list_atom_surf[*]}"
+
+echo ' -- ADSORPTION SUTRUCTURES HAVE BEEN CALCULATED --'
+
+# Changes labels of elements to its original definition, thus allowing different kinds of same element (Fe1 Fe2) with different properties (e.g. spin state in compounds with magnetic coupling) 
+
+${DockOnSurf_path}/modules/fe_change.sh $molecule $surface ${MolOnSurf_results_path}/${molecule}
+
+### Lancement des calculs pour les molecules adsorbees
+
+${DockOnSurf_path}/modules/launch_script+diss.sh $molecule
+
+### Attente jusqu'à ce que tous les calculs aient fini
+
+mol2="$(echo $molecule | cut -c1-5)"
+while [ "$(qstat | grep surf_${mol2} | wc -l)" != 0 ];
+do 
+	sleep 300
+done
+
+echo ' -- ALL CALCULATIONS FOR ADSORPTION STRUCTURES HAVE BEEN DONE --'
+
+### Suppression des fichiers inutiles
+
+rm -r ${MolOnSurf_results_path}/${molecule}/*/*RESTART*
+rm -r ${MolOnSurf_results_path}/${molecule}/*/*BFGS*
+rm -r ${MolOnSurf_results_path}/${molecule}/*/*.p*
+
+### Clustering et lancement des structures centres des clusters
+
+nb_at_surf="$(awk 'NR==1{print}' ${Surface_path}/${surface}.xyz)"
+num_center="$(expr ${nb_at_surf} + $center )"
+
+${DockOnSurf_path}/modules/script_grandes_molecules+diss.sh $molecule $num_center ${nb_at_surf} ${cutoff}
+
+echo ' -- CLUSTERING HAS BEEN DONE AND CLUSTER CENTERS CALCULATIONS HAVE BEEN LAUNCHED --'
+
+### Attente jusqu'à ce que tous les calculs aient fini
+
+while [ "$(qstat | grep surf_${mol2} | wc -l)" != 0 ];
+do
+        sleep 300
+done
+
+echo ' -- ALL CLUSTER CENTERS CALCULATIONS HAVE BEEN DONE --'
+
+
+echo " "
+echo ' -- RESULTS --'
+
+echo ' '
+echo "Most stable structure for molecule alone: ${numero_mol_stable} with energy ${min_energy_molecule} Ha"
+
+min_energy_global=0
+for global_struct in ${MolOnSurf_results_path}/${molecule}/relaunched_calculations/* ; do
+	if [[ "$global_struct" != *.inp ]]
+	then
+		if [[ "$global_struct" != *.j ]]
+		then 
+			energy_global="$(awk '/ENERGY/{E=$NF} END{print E}' ${global_struct}/*${molecule}.out)"
+			if [[ "${energy_global}" > "${min_energy_global}" ]]
+			then
+				min_energy_global=${energy_global}
+				numero_global_stable="$(echo ${global_struct} | awk -F'/' '{print $NF}')"
+			fi
+		fi
+	fi
+done
+
+echo "Most stable structure for adsorption: ${numero_global_stable} with energy ${min_energy_global} Ha"
+
+
diff --git a/modules/add_adsorbate_euler+dissociation.py b/modules/add_adsorbate_euler+dissociation.py
new file mode 100755
index 0000000..3ea059a
--- /dev/null
+++ b/modules/add_adsorbate_euler+dissociation.py
@@ -0,0 +1,278 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+import ase
+from ase import io
+from ase.io import read
+from ase.io import write
+from ase import Atom
+from ase import build
+from ase.build import add_adsorbate
+from ase import constraints
+from ase.constraints import FixAtoms
+from ase import neighborlist
+from ase.neighborlist import NeighborList
+from ase import utils
+import ase_Fe1_Fe2
+
+if len(sys.argv) < 7 :
+    print "This script has to be used with at least 6 arguments."
+    print "1. Name of the file containing the surface"
+    print "2. Atom number of the atom of the surface where the adsorbate will be adsorbed"
+    print "3. Name of the file containing the molecule to be adsorbed"
+    print "4. Atom number of the atom of the molecule that will be adsorbed of the surface"
+    print "5. Distance between the molecule and the surface"
+    print "6. Name of the output file"
+    print "Other arguments can be added before the name of the output file :"
+    print "  Second atom of the molecule"
+    print "  phi angle (rotation around z)"
+    print "  theta angle (rotation around new x)"
+    print "  psi angle (rotation around new z) "
+    sys.exit(1)
+
+
+#############################################################
+# Définition de la fonction permettant de calculer les angles
+#############################################################
+
+def get_proper_angle(v1, v2, degrees=True):
+    norm_dot = np.dot(v1/np.linalg.norm(v1), v2/np.linalg.norm(v2))
+    angle = np.arccos(np.sign(norm_dot) if np.abs(norm_dot) >= 1 else norm_dot)
+    return(angle*180/np.pi if degrees else angle)
+
+
+##########################################
+# Lecture des différents arguments de base
+##########################################
+
+surface_file = sys.argv[1]
+atom_surface = int(sys.argv[2])
+molecule_file = sys.argv[3]
+atom_molecule = int(sys.argv[4])
+distance = int(sys.argv[5])
+if len(sys.argv) < 8 :
+    output = sys.argv[6]
+
+surface = read(surface_file)
+molecule = read(molecule_file)
+
+#################################
+# Lecture des positions de atomes
+#################################
+
+atom_surf_1 = surface[atom_surface].position
+x_atom_surf_1 = atom_surf_1[0]
+y_atom_surf_1 = atom_surf_1[1]
+z_atom_surf_1 = atom_surf_1[2]
+
+atom_mol_1 = molecule[atom_molecule].position
+x_atom_mol_1 = atom_mol_1[0]
+y_atom_mol_1 = atom_mol_1[1]
+z_atom_mol_1 = atom_mol_1[2]
+
+nb_atom_mol = int(len(molecule))
+
+
+#######################################################
+# Cas avec un plus grand nombre d'arguments (rotations)
+#######################################################
+
+
+if len(sys.argv) > 7 :
+
+    #######################
+    # Lecture des arguments
+    #######################
+
+    atom_molecule_2 = int(sys.argv[6])
+    phi = int(sys.argv[7])
+    theta = int(sys.argv[8])
+    psi = int(sys.argv[9])
+    output = sys.argv[10]
+
+    theta_min = theta - 0.01
+    theta_max = theta + 0.01
+    phi_min = phi - 0.01
+    phi_max = phi + 0.01
+    psi_min = psi - 0.01
+    psi_max = psi + 0.01
+
+    atom_mol_2 = molecule[atom_molecule_2].position
+    x_atom_mol_2 = atom_mol_2[0]
+    y_atom_mol_2 = atom_mol_2[1]
+    z_atom_mol_2 = atom_mol_2[2]
+
+    ##################################
+    # Rotation avec les angles d'Euler
+    ##################################
+
+    #Initialisation : placement du second atome de la molecule sur Oy
+    vect_y = [0, 1, 0]
+    vect_atom_mol_2 = molecule[atom_molecule_2].position - molecule[atom_molecule].position
+    angle_1 = get_proper_angle(vect_y, vect_atom_mol_2)
+    if (angle_1 != 0) :
+        vect_normal = np.cross(vect_y, vect_atom_mol_2)
+        molecule.rotate(-angle_1, vect_normal, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
+        vect_atom_mol_2_verif = molecule[atom_molecule_2].position - molecule[atom_molecule].position
+        angle_1_verif = get_proper_angle(vect_atom_mol_2_verif, vect_y)
+        angle_max = 0.01
+        angle_min = -0.01
+        if (angle_1_verif < angle_min or angle_1_verif > angle_max) :
+            print 'Error in initialisation'
+
+    #Rotation Euler
+    molecule.euler_rotate(phi, theta, psi, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
+
+    #Correction des collisions entre la molécule et la surface
+    z_atom_surf_max = z_atom_surf_1 + 1
+
+    collision = 0.5
+    rotation_tot = 0
+    nb_atom_collision_min = 0
+    rotation_opt = 0
+    rotation_modif = 5
+    
+    for i in range(0,nb_atom_mol) :
+        atom_test = molecule[i].position
+        z_atom_test = atom_test[2]
+        z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
+        if (z_atom_test_final < z_atom_surf_max) :
+            collision = 1
+            nb_atom_collision_min += 1
+    if collision == 1 :
+        print 'Collision between the molecule and the surface - modification of theta angle'
+
+    vect_z = [0, 0, 1]
+
+    while (collision == 1 and rotation_tot < 354) :
+        vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
+        
+        if (vect_atom_mol_2_post_euler[0] != 0 and vect_atom_mol_2_post_euler[1] != 0) :
+            vect_rotation_theta = np.cross(vect_atom_mol_2_post_euler, vect_z)
+        else :
+            vect_rotation_theta = [0.0001, 0.0001 , 0]
+        molecule.rotate(rotation_modif, vect_rotation_theta, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
+        rotation_tot += 5
+        
+        collision = 0
+        nb_atom_collision = 0
+
+        for i in range(0,nb_atom_mol) :
+            atom_test = molecule[i].position
+            z_atom_test = atom_test[2]
+            z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
+            if (z_atom_test_final < z_atom_surf_max) :
+                nb_atom_collision += 1
+                collision = 1
+        if (nb_atom_collision < nb_atom_collision_min) :
+            nb_atom_collision_min = nb_atom_collision
+            rotation_opt = rotation_tot
+
+        if nb_atom_collision_min != 0 :
+            rotation_tot_2 = 0
+            rotation_opt_2 = 0
+
+            while (collision == 1 and rotation_tot_2 < 354) :
+                vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
+                molecule.rotate(rotation_modif, vect_atom_mol_2_post_euler, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
+                rotation_tot_2 += 5
+                
+                collision = 0
+                nb_atom_collision = 0
+
+                for i in range(0,nb_atom_mol) :
+                    atom_test = molecule[i].position
+                    z_atom_test = atom_test[2]
+                    z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
+                    if (z_atom_test_final < z_atom_surf_max) :
+                        nb_atom_collision += 1
+                        collision = 1
+                if (nb_atom_collision < nb_atom_collision_min) :
+                    nb_atom_collision_min = nb_atom_collision
+                    rotation_opt_2 = rotation_tot_2
+
+    if collision == 0 :
+        print 'Collision corrected'
+    elif collision == 1 :
+        print 'Error: the collision could not be corrected'
+        pos1 = output.rfind('/')
+        pos = pos1 + 1
+        num_coll=output[pos:]
+        j = output[:pos1]
+        pos2 = j.rfind('/')
+        molecule_directory = j[:pos2]
+        fichier = open("%s/errors" % molecule_directory , "a")
+        fichier.write("Adsorption %s : ERROR the collision could not be corrected\n" % num_coll )
+        fichier.close()
+
+
+##########################################
+# Adsorption de la molecule sur la surface
+##########################################
+
+add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1), mol_index=atom_molecule)
+
+out=output+".xyz"
+write(out, surface)
+
+
+############################
+# Dissociation si necessaire
+############################
+
+if collision == 0 :
+    surface = read(surface_file)
+    molecule_reminder = molecule
+    
+    list_cutoffs = utils.natural_cutoffs(molecule)
+    nl=NeighborList(list_cutoffs, self_interaction=False, bothways=True)
+
+    nl.update(molecule)
+    neighbor_indices, offsets = nl.get_neighbors(atom_molecule)
+
+    symbols = molecule.get_chemical_symbols()
+
+    neighbors_symbols=[]
+    for i in neighbor_indices:
+            neighbors_symbols.append(symbols[i])
+
+    nb_neighbors = len(neighbor_indices)
+
+    diss=0
+    for i in range(0,nb_neighbors):
+        if neighbors_symbols[i] == "H":
+            diss=1
+            H_atom=neighbor_indices[i] #nb of the H atom
+            print "Dissociation possible"
+
+    if diss == 1:
+        list_cutoffs_surface = utils.natural_cutoffs(surface)
+        nl2=NeighborList(list_cutoffs_surface, self_interaction=False, bothways=True)
+        nl2.update(surface)
+        neighbor_indices_surf, offsets_surf = nl2.get_neighbors(atom_surface)
+        symbols_surf = surface.get_chemical_symbols()
+        for i in neighbor_indices_surf:
+            surface = read(surface_file)
+            molecule = molecule_reminder
+            neighbor_surf_coord = surface[i].position
+            coord_H_atom=molecule[H_atom].position
+            if (neighbor_surf_coord[2] > z_atom_surf_1-1 and symbols_surf[i] == "O"):
+                vector_H_diss=[neighbor_surf_coord[0]-coord_H_atom[0], neighbor_surf_coord[1]-coord_H_atom[1], neighbor_surf_coord[2]+1-coord_H_atom[2]]
+                translation_matrix_dissociation=[]
+                for k in range (0, nb_atom_mol):
+                    if k != H_atom:
+                        translation_matrix_dissociation.append([0,0,0])
+                    else:
+                        translation_matrix_dissociation.append(vector_H_diss)
+                molecule.translate(translation_matrix_dissociation)
+                add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1), mol_index=atom_molecule)
+                output_diss = output + "_diss_" + str(i) + ".xyz"
+                write(output_diss, surface)
+    
+
+
+
+
+
diff --git a/modules/ase_Fe1_Fe2.py b/modules/ase_Fe1_Fe2.py
new file mode 100755
index 0000000..46090a5
--- /dev/null
+++ b/modules/ase_Fe1_Fe2.py
@@ -0,0 +1,58 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+
+import ase
+from ase import Atom
+from ase import data
+import numpy as np
+
+# Chemical symbols of new types of atoms
+data.chemical_symbols += ['Fe1', 'Fe2']
+data.chemical_symbols += ['O1']
+
+# Atomic numbers
+Z_Fe = data.atomic_numbers['Fe']
+Z_O = data.atomic_numbers['O']
+data.atomic_numbers['Fe1'] = Z_Fe
+data.atomic_numbers['Fe2'] = Z_Fe
+data.atomic_numbers['O1'] = Z_O
+
+# Atomic names
+data.atomic_names += ['Iron', 'Iron']
+data.atomic_names += ['Oxygen'] 
+
+# Atomic masses
+Fe_mass = data.atomic_masses_iupac2016[Z_Fe]
+O_mass = data.atomic_masses_iupac2016[Z_O] 
+data.atomic_masses_iupac2016 = np.concatenate((data.atomic_masses_iupac2016, np.array([Fe_mass,Fe_mass])))
+data.atomic_masses_iupac2016 = np.concatenate((data.atomic_masses_iupac2016, np.array([O_mass])))
+
+data.atomic_masses = data.atomic_masses_iupac2016
+
+Fe_mass = data.atomic_masses_legacy[Z_Fe]
+O_mass = data.atomic_masses_legacy[Z_O]
+data.atomic_masses_legacy = np.concatenate((data.atomic_masses_legacy, np.array([Fe_mass,Fe_mass])))
+data.atomic_masses_legacy = np.concatenate((data.atomic_masses_legacy, np.array([O_mass])))
+
+# Atomic radius
+Fe_radius = data.covalent_radii[Z_Fe]
+O_radius = data.covalent_radii[Z_O]
+data.covalent_radii = np.concatenate((data.covalent_radii, np.array([Fe_radius,Fe_radius])))
+data.covalent_radii = np.concatenate((data.covalent_radii, np.array([O_radius])))
+
+# Reference and ground state datas
+data.reference_states += [data.reference_states[Z_Fe], data.reference_states[Z_Fe]]
+data.reference_states += [data.reference_states[Z_O]]
+
+Fe_gsmm = data.ground_state_magnetic_moments[Z_Fe]
+O_gsmm = data.ground_state_magnetic_moments[Z_O]
+
+data.ground_state_magnetic_moments = np.concatenate((data.ground_state_magnetic_moments, np.array([Fe_gsmm,Fe_gsmm])))
+data.ground_state_magnetic_moments = np.concatenate((data.ground_state_magnetic_moments, np.array([O_gsmm])))
+
+# VdW raduis
+Fe_vdw_r = data.vdw_radii[Z_Fe]
+O_vdw_r = data.vdw_radii[Z_O]
+data.vdw_radii = np.concatenate((data.vdw_radii, np.array([Fe_vdw_r,Fe_vdw_r])))
+data.vdw_radii = np.concatenate((data.vdw_radii, np.array([O_vdw_r])))
diff --git a/modules/calculate_rmsd b/modules/calculate_rmsd
new file mode 100755
index 0000000..65c9212
--- /dev/null
+++ b/modules/calculate_rmsd
@@ -0,0 +1,1042 @@
+#!/usr/bin/env python
+__doc__ = \
+"""
+Calculate Root-mean-square deviation (RMSD) between structure A and B, in XYZ
+or PDB format, using transformation and rotation.
+
+For more information, usage, example and citation read more at
+https://github.com/charnley/rmsd
+"""
+
+__version__ = '1.3.2'
+
+import copy
+import re
+
+import numpy as np
+from scipy.optimize import linear_sum_assignment
+from scipy.spatial.distance import cdist
+
+
+AXIS_SWAPS = np.array([
+    [0, 1, 2],
+    [0, 2, 1],
+    [1, 0, 2],
+    [1, 2, 0],
+    [2, 1, 0],
+    [2, 0, 1]])
+
+AXIS_REFLECTIONS = np.array([
+    [1, 1, 1],
+    [-1, 1, 1],
+    [1, -1, 1],
+    [1, 1, -1],
+    [-1, -1, 1],
+    [-1, 1, -1],
+    [1, -1, -1],
+    [-1, -1, -1]])
+
+
+def rmsd(V, W):
+    """
+    Calculate Root-mean-square deviation from two sets of vectors V and W.
+
+    Parameters
+    ----------
+    V : array
+        (N,D) matrix, where N is points and D is dimension.
+    W : array
+        (N,D) matrix, where N is points and D is dimension.
+
+    Returns
+    -------
+    rmsd : float
+        Root-mean-square deviation between the two vectors
+    """
+    D = len(V[0])
+    N = len(V)
+    result = 0.0
+    for v, w in zip(V, W):
+        result += sum([(v[i] - w[i])**2.0 for i in range(D)])
+    return np.sqrt(result/N)
+
+
+def kabsch_rmsd(P, Q, translate=False):
+    """
+    Rotate matrix P unto Q using Kabsch algorithm and calculate the RMSD.
+
+    Parameters
+    ----------
+    P : array
+        (N,D) matrix, where N is points and D is dimension.
+    Q : array
+        (N,D) matrix, where N is points and D is dimension.
+    translate : bool
+        Use centroids to translate vector P and Q unto each other.
+
+    Returns
+    -------
+    rmsd : float
+        root-mean squared deviation
+    """
+    if translate:
+        Q = Q - centroid(Q)
+        P = P - centroid(P)
+
+    P = kabsch_rotate(P, Q)
+    return rmsd(P, Q)
+
+
+def kabsch_rotate(P, Q):
+    """
+    Rotate matrix P unto matrix Q using Kabsch algorithm.
+
+    Parameters
+    ----------
+    P : array
+        (N,D) matrix, where N is points and D is dimension.
+    Q : array
+        (N,D) matrix, where N is points and D is dimension.
+
+    Returns
+    -------
+    P : array
+        (N,D) matrix, where N is points and D is dimension,
+        rotated
+
+    """
+    U = kabsch(P, Q)
+
+    # Rotate P
+    P = np.dot(P, U)
+    return P
+
+
+def kabsch(P, Q):
+    """
+    Using the Kabsch algorithm with two sets of paired point P and Q, centered
+    around the centroid. Each vector set is represented as an NxD
+    matrix, where D is the the dimension of the space.
+
+    The algorithm works in three steps:
+    - a centroid translation of P and Q (assumed done before this function
+      call)
+    - the computation of a covariance matrix C
+    - computation of the optimal rotation matrix U
+
+    For more info see http://en.wikipedia.org/wiki/Kabsch_algorithm
+
+    Parameters
+    ----------
+    P : array
+        (N,D) matrix, where N is points and D is dimension.
+    Q : array
+        (N,D) matrix, where N is points and D is dimension.
+
+    Returns
+    -------
+    U : matrix
+        Rotation matrix (D,D)
+    """
+
+    # Computation of the covariance matrix
+    C = np.dot(np.transpose(P), Q)
+
+    # Computation of the optimal rotation matrix
+    # This can be done using singular value decomposition (SVD)
+    # Getting the sign of the det(V)*(W) to decide
+    # whether we need to correct our rotation matrix to ensure a
+    # right-handed coordinate system.
+    # And finally calculating the optimal rotation matrix U
+    # see http://en.wikipedia.org/wiki/Kabsch_algorithm
+    V, S, W = np.linalg.svd(C)
+    d = (np.linalg.det(V) * np.linalg.det(W)) < 0.0
+
+    if d:
+        S[-1] = -S[-1]
+        V[:, -1] = -V[:, -1]
+
+    # Create Rotation matrix U
+    U = np.dot(V, W)
+
+    return U
+
+
+def quaternion_rmsd(P, Q):
+    """
+    Rotate matrix P unto Q and calculate the RMSD
+    based on doi:10.1016/1049-9660(91)90036-O
+
+    Parameters
+    ----------
+    P : array
+        (N,D) matrix, where N is points and D is dimension.
+    Q : array
+        (N,D) matrix, where N is points and D is dimension.
+
+    Returns
+    -------
+    rmsd : float
+    """
+    rot = quaternion_rotate(P, Q)
+    P = np.dot(P, rot)
+    return rmsd(P, Q)
+
+
+def quaternion_transform(r):
+    """
+    Get optimal rotation
+    note: translation will be zero when the centroids of each molecule are the
+    same
+    """
+    Wt_r = makeW(*r).T
+    Q_r = makeQ(*r)
+    rot = Wt_r.dot(Q_r)[:3, :3]
+    return rot
+
+
+def makeW(r1, r2, r3, r4=0):
+    """
+    matrix involved in quaternion rotation
+    """
+    W = np.asarray([
+        [r4, r3, -r2, r1],
+        [-r3, r4, r1, r2],
+        [r2, -r1, r4, r3],
+        [-r1, -r2, -r3, r4]])
+    return W
+
+
+def makeQ(r1, r2, r3, r4=0):
+    """
+    matrix involved in quaternion rotation
+    """
+    Q = np.asarray([
+        [r4, -r3, r2, r1],
+        [r3, r4, -r1, r2],
+        [-r2, r1, r4, r3],
+        [-r1, -r2, -r3, r4]])
+    return Q
+
+
+def quaternion_rotate(X, Y):
+    """
+    Calculate the rotation
+
+    Parameters
+    ----------
+    X : array
+        (N,D) matrix, where N is points and D is dimension.
+    Y: array
+        (N,D) matrix, where N is points and D is dimension.
+
+    Returns
+    -------
+    rot : matrix
+        Rotation matrix (D,D)
+    """
+    N = X.shape[0]
+    W = np.asarray([makeW(*Y[k]) for k in range(N)])
+    Q = np.asarray([makeQ(*X[k]) for k in range(N)])
+    Qt_dot_W = np.asarray([np.dot(Q[k].T, W[k]) for k in range(N)])
+    W_minus_Q = np.asarray([W[k] - Q[k] for k in range(N)])
+    A = np.sum(Qt_dot_W, axis=0)
+    eigen = np.linalg.eigh(A)
+    r = eigen[1][:, eigen[0].argmax()]
+    rot = quaternion_transform(r)
+    return rot
+
+
+def centroid(X):
+    """
+    Centroid is the mean position of all the points in all of the coordinate
+    directions, from a vectorset X.
+
+    https://en.wikipedia.org/wiki/Centroid
+
+    C = sum(X)/len(X)
+
+    Parameters
+    ----------
+    X : array
+        (N,D) matrix, where N is points and D is dimension.
+
+    Returns
+    -------
+    C : float
+        centroid
+    """
+    C = X.mean(axis=0)
+    return C
+
+
+def reorder_distance(p_atoms, q_atoms, p_coord, q_coord):
+    """
+    Re-orders the input atom list and xyz coordinates by atom type and then by
+    distance of each atom from the centroid.
+
+    Parameters
+    ----------
+    atoms : array
+        (N,1) matrix, where N is points holding the atoms' names
+    coord : array
+        (N,D) matrix, where N is points and D is dimension
+
+    Returns
+    -------
+    atoms_reordered : array
+        (N,1) matrix, where N is points holding the ordered atoms' names
+    coords_reordered : array
+        (N,D) matrix, where N is points and D is dimension (rows re-ordered)
+    """
+
+    # Find unique atoms
+    unique_atoms = np.unique(p_atoms)
+
+    # generate full view from q shape to fill in atom view on the fly
+    view_reorder = np.zeros(q_atoms.shape, dtype=int)
+
+    for atom in unique_atoms:
+
+        p_atom_idx, = np.where(p_atoms == atom)
+        q_atom_idx, = np.where(q_atoms == atom)
+
+        A_coord = p_coord[p_atom_idx]
+        B_coord = q_coord[q_atom_idx]
+
+        # Calculate distance from each atom to centroid
+        A_norms = np.linalg.norm(A_coord, axis=1)
+        B_norms = np.linalg.norm(B_coord, axis=1)
+
+        reorder_indices_A = np.argsort(A_norms)
+        reorder_indices_B = np.argsort(B_norms)
+
+        # Project the order of P onto Q
+        translator = np.argsort(reorder_indices_A)
+        view = reorder_indices_B[translator]
+        view_reorder[p_atom_idx] = q_atom_idx[view]
+
+    return view_reorder
+
+
+def hungarian(A, B):
+    """
+    Hungarian reordering.
+
+    Assume A and B are coordinates for atoms of SAME type only
+    """
+
+    # should be kabasch here i think
+    distances = cdist(A, B, 'euclidean')
+
+    # Perform Hungarian analysis on distance matrix between atoms of 1st
+    # structure and trial structure
+    indices_a, indices_b = linear_sum_assignment(distances)
+
+    return indices_b
+
+
+def reorder_hungarian(p_atoms, q_atoms, p_coord, q_coord):
+    """
+    Re-orders the input atom list and xyz coordinates using the Hungarian
+    method (using optimized column results)
+
+    Parameters
+    ----------
+    p_atoms : array
+        (N,1) matrix, where N is points holding the atoms' names
+    p_atoms : array
+        (N,1) matrix, where N is points holding the atoms' names
+    p_coord : array
+        (N,D) matrix, where N is points and D is dimension
+    q_coord : array
+        (N,D) matrix, where N is points and D is dimension
+
+    Returns
+    -------
+    view_reorder : array
+             (N,1) matrix, reordered indexes of atom alignment based on the
+             coordinates of the atoms
+
+    """
+
+    # Find unique atoms
+    unique_atoms = np.unique(p_atoms)
+
+    # generate full view from q shape to fill in atom view on the fly
+    view_reorder = np.zeros(q_atoms.shape, dtype=int)
+    view_reorder -= 1
+
+    for atom in unique_atoms:
+        p_atom_idx, = np.where(p_atoms == atom)
+        q_atom_idx, = np.where(q_atoms == atom)
+
+        A_coord = p_coord[p_atom_idx]
+        B_coord = q_coord[q_atom_idx]
+
+        view = hungarian(A_coord, B_coord)
+        view_reorder[p_atom_idx] = q_atom_idx[view]
+
+    return view_reorder
+
+
+def generate_permutations(elements, n):
+    """
+    Heap's algorithm for generating all n! permutations in a list
+    https://en.wikipedia.org/wiki/Heap%27s_algorithm
+
+    """
+    c = [0] * n
+    yield elements
+    i = 0
+    while i < n:
+        if c[i] < i:
+            if i % 2 == 0:
+                elements[0], elements[i] = elements[i], elements[0]
+            else:
+                elements[c[i]], elements[i] = elements[i], elements[c[i]]
+            yield elements
+            c[i] += 1
+            i = 0
+        else:
+            c[i] = 0
+            i += 1
+
+
+def brute_permutation(A, B):
+    """
+    Re-orders the input atom list and xyz coordinates using the brute force
+    method of permuting all rows of the input coordinates
+
+    Parameters
+    ----------
+    A : array
+        (N,D) matrix, where N is points and D is dimension
+    B : array
+        (N,D) matrix, where N is points and D is dimension
+
+    Returns
+    -------
+    view : array
+        (N,1) matrix, reordered view of B projected to A
+    """
+
+    rmsd_min = np.inf
+    view_min = None
+
+    # Sets initial ordering for row indices to [0, 1, 2, ..., len(A)], used in
+    # brute-force method
+
+    num_atoms = A.shape[0]
+    initial_order = list(range(num_atoms))
+
+    for reorder_indices in generate_permutations(initial_order, num_atoms):
+
+        # Re-order the atom array and coordinate matrix
+        coords_ordered = B[reorder_indices]
+
+        # Calculate the RMSD between structure 1 and the Hungarian re-ordered
+        # structure 2
+        rmsd_temp = kabsch_rmsd(A, coords_ordered)
+
+        # Replaces the atoms and coordinates with the current structure if the
+        # RMSD is lower
+        if rmsd_temp < rmsd_min:
+            rmsd_min = rmsd_temp
+            view_min = copy.deepcopy(reorder_indices)
+
+    return view_min
+
+
+def reorder_brute(p_atoms, q_atoms, p_coord, q_coord):
+    """
+    Re-orders the input atom list and xyz coordinates using all permutation of
+    rows (using optimized column results)
+
+    Parameters
+    ----------
+    p_atoms : array
+        (N,1) matrix, where N is points holding the atoms' names
+    q_atoms : array
+        (N,1) matrix, where N is points holding the atoms' names
+    p_coord : array
+        (N,D) matrix, where N is points and D is dimension
+    q_coord : array
+        (N,D) matrix, where N is points and D is dimension
+
+    Returns
+    -------
+    view_reorder : array
+        (N,1) matrix, reordered indexes of atom alignment based on the
+        coordinates of the atoms
+
+    """
+
+    # Find unique atoms
+    unique_atoms = np.unique(p_atoms)
+
+    # generate full view from q shape to fill in atom view on the fly
+    view_reorder = np.zeros(q_atoms.shape, dtype=int)
+    view_reorder -= 1
+
+    for atom in unique_atoms:
+        p_atom_idx, = np.where(p_atoms == atom)
+        q_atom_idx, = np.where(q_atoms == atom)
+
+        A_coord = p_coord[p_atom_idx]
+        B_coord = q_coord[q_atom_idx]
+
+        view = brute_permutation(A_coord, B_coord)
+        view_reorder[p_atom_idx] = q_atom_idx[view]
+
+    return view_reorder
+
+
+def check_reflections(p_atoms, q_atoms, p_coord, q_coord,
+                      reorder_method=reorder_hungarian,
+                      rotation_method=kabsch_rmsd,
+                      keep_stereo=False):
+    """
+    Minimize RMSD using reflection planes for molecule P and Q
+
+    Warning: This will affect stereo-chemistry
+
+    Parameters
+    ----------
+    p_atoms : array
+        (N,1) matrix, where N is points holding the atoms' names
+    q_atoms : array
+        (N,1) matrix, where N is points holding the atoms' names
+    p_coord : array
+        (N,D) matrix, where N is points and D is dimension
+    q_coord : array
+        (N,D) matrix, where N is points and D is dimension
+
+    Returns
+    -------
+    min_rmsd
+    min_swap
+    min_reflection
+    min_review
+
+    """
+
+    min_rmsd = np.inf
+    min_swap = None
+    min_reflection = None
+    min_review = None
+    tmp_review = None
+    swap_mask = [1,-1,-1,1,-1,1]
+    reflection_mask = [1,-1,-1,-1,1,1,1,-1]
+
+    for swap, i in zip(AXIS_SWAPS, swap_mask):
+        for reflection, j in zip(AXIS_REFLECTIONS, reflection_mask):
+            if keep_stereo and  i * j == -1: continue # skip enantiomers
+
+            tmp_atoms = copy.copy(q_atoms)
+            tmp_coord = copy.deepcopy(q_coord)
+            tmp_coord = tmp_coord[:, swap]
+            tmp_coord = np.dot(tmp_coord, np.diag(reflection))
+            tmp_coord -= centroid(tmp_coord)
+
+            # Reorder
+            if reorder_method is not None:
+                tmp_review = reorder_method(p_atoms, tmp_atoms, p_coord, tmp_coord)
+                tmp_coord = tmp_coord[tmp_review]
+                tmp_atoms = tmp_atoms[tmp_review]
+
+            # Rotation
+            if rotation_method is None:
+                this_rmsd = rmsd(p_coord, tmp_coord)
+            else:
+                this_rmsd = rotation_method(p_coord, tmp_coord)
+
+            if this_rmsd < min_rmsd:
+                min_rmsd = this_rmsd
+                min_swap = swap
+                min_reflection = reflection
+                min_review = tmp_review
+
+    if not (p_atoms == q_atoms[min_review]).all():
+        print("error: Not aligned")
+        quit()
+
+    return min_rmsd, min_swap, min_reflection, min_review
+
+
+def set_coordinates(atoms, V, title="", decimals=8):
+    """
+    Print coordinates V with corresponding atoms to stdout in XYZ format.
+    Parameters
+    ----------
+    atoms : list
+        List of atomic types
+    V : array
+        (N,3) matrix of atomic coordinates
+    title : string (optional)
+        Title of molecule
+    decimals : int (optional)
+        number of decimals for the coordinates
+
+    Return
+    ------
+    output : str
+        Molecule in XYZ format
+
+    """
+    N, D = V.shape
+
+    fmt = "{:2s}" + (" {:15."+str(decimals)+"f}")*3
+
+    out = list()
+    out += [str(N)]
+    out += [title]
+
+    for i in range(N):
+        atom = atoms[i]
+        atom = atom[0].upper() + atom[1:]
+        out += [fmt.format(atom, V[i, 0], V[i, 1], V[i, 2])]
+
+    return "\n".join(out)
+
+
+def print_coordinates(atoms, V, title=""):
+    """
+    Print coordinates V with corresponding atoms to stdout in XYZ format.
+
+    Parameters
+    ----------
+    atoms : list
+        List of element types
+    V : array
+        (N,3) matrix of atomic coordinates
+    title : string (optional)
+        Title of molecule
+
+    """
+
+    print(set_coordinates(atoms, V, title=title))
+
+    return
+
+
+def get_coordinates(filename, fmt):
+    """
+    Get coordinates from filename in format fmt. Supports XYZ and PDB.
+    Parameters
+    ----------
+    filename : string
+        Filename to read
+    fmt : string
+        Format of filename. Either xyz or pdb.
+    Returns
+    -------
+    atoms : list
+        List of atomic types
+    V : array
+        (N,3) where N is number of atoms
+    """
+    if fmt == "xyz":
+        get_func = get_coordinates_xyz
+    elif fmt == "pdb":
+        get_func = get_coordinates_pdb
+    else:
+        exit("Could not recognize file format: {:s}".format(fmt))
+
+    return get_func(filename)
+
+
+def get_coordinates_pdb(filename):
+    """
+    Get coordinates from the first chain in a pdb file
+    and return a vectorset with all the coordinates.
+
+    Parameters
+    ----------
+    filename : string
+        Filename to read
+
+    Returns
+    -------
+    atoms : list
+        List of atomic types
+    V : array
+        (N,3) where N is number of atoms
+    """
+
+    # PDB files tend to be a bit of a mess. The x, y and z coordinates
+    # are supposed to be in column 31-38, 39-46 and 47-54, but this is
+    # not always the case.
+    # Because of this the three first columns containing a decimal is used.
+    # Since the format doesn't require a space between columns, we use the
+    # above column indices as a fallback.
+
+    x_column = None
+    V = list()
+
+    # Same with atoms and atom naming.
+    # The most robust way to do this is probably
+    # to assume that the atomtype is given in column 3.
+
+    atoms = list()
+
+    with open(filename, 'r') as f:
+        lines = f.readlines()
+        for line in lines:
+            if line.startswith("TER") or line.startswith("END"):
+                break
+            if line.startswith("ATOM"):
+                tokens = line.split()
+                # Try to get the atomtype
+                try:
+                    atom = tokens[2][0]
+                    if atom in ("H", "C", "N", "O", "S", "P"):
+                        atoms.append(atom)
+                    else:
+                        # e.g. 1HD1
+                        atom = tokens[2][1]
+                        if atom == "H":
+                            atoms.append(atom)
+                        else:
+                            raise Exception
+                except:
+                    exit("error: Parsing atomtype for the following line: \n{0:s}".format(line))
+
+                if x_column == None:
+                    try:
+                        # look for x column
+                        for i, x in enumerate(tokens):
+                            if "." in x and "." in tokens[i + 1] and "." in tokens[i + 2]:
+                                x_column = i
+                                break
+                    except IndexError:
+                        exit("error: Parsing coordinates for the following line: \n{0:s}".format(line))
+                # Try to read the coordinates
+                try:
+                    V.append(np.asarray(tokens[x_column:x_column + 3], dtype=float))
+                except:
+                    # If that doesn't work, use hardcoded indices
+                    try:
+                        x = line[30:38]
+                        y = line[38:46]
+                        z = line[46:54]
+                        V.append(np.asarray([x, y ,z], dtype=float))
+                    except:
+                        exit("error: Parsing input for the following line: \n{0:s}".format(line))
+
+
+    V = np.asarray(V)
+    atoms = np.asarray(atoms)
+
+    assert V.shape[0] == atoms.size
+
+    return atoms, V
+
+
+def get_coordinates_xyz(filename):
+    """
+    Get coordinates from filename and return a vectorset with all the
+    coordinates, in XYZ format.
+
+    Parameters
+    ----------
+    filename : string
+        Filename to read
+
+    Returns
+    -------
+    atoms : list
+        List of atomic types
+    V : array
+        (N,3) where N is number of atoms
+    """
+
+    f = open(filename, 'r')
+    V = list()
+    atoms = list()
+    n_atoms = 0
+
+    # Read the first line to obtain the number of atoms to read
+    try:
+        n_atoms = int(f.readline())
+    except ValueError:
+        exit("error: Could not obtain the number of atoms in the .xyz file.")
+
+    # Skip the title line
+    f.readline()
+
+    # Use the number of atoms to not read beyond the end of a file
+    for lines_read, line in enumerate(f):
+
+        if lines_read == n_atoms:
+            break
+
+        atom = re.findall(r'[a-zA-Z]+', line)[0]
+        atom = atom.upper()
+
+        numbers = re.findall(r'[-]?\d+\.\d*(?:[Ee][-\+]\d+)?', line)
+        numbers = [float(number) for number in numbers]
+
+        # The numbers are not valid unless we obtain exacly three
+        if len(numbers) >= 3:
+            V.append(np.array(numbers)[:3])
+            atoms.append(atom)
+        else:
+            exit("Reading the .xyz file failed in line {0}. Please check the format.".format(lines_read + 2))
+
+    f.close()
+    atoms = np.array(atoms)
+    V = np.array(V)
+    return atoms, V
+
+
+def main():
+
+    import argparse
+    import sys
+
+    description = __doc__
+
+    version_msg = """
+rmsd {}
+
+See https://github.com/charnley/rmsd for citation information
+
+"""
+    version_msg = version_msg.format(__version__)
+
+    epilog = """
+"""
+
+    parser = argparse.ArgumentParser(
+        usage='calculate_rmsd [options] FILE_A FILE_B',
+        description=description,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog=epilog)
+
+
+    # Input structures
+    parser.add_argument('structure_a', metavar='FILE_A', type=str, help='structures in .xyz or .pdb format')
+    parser.add_argument('structure_b', metavar='FILE_B', type=str)
+
+    # Admin
+    parser.add_argument('-v', '--version', action='version', version=version_msg)
+
+    # Rotation
+    parser.add_argument('-r', '--rotation', action='store', default="kabsch", help='select rotation method. "kabsch" (default), "quaternion" or "none"', metavar="METHOD")
+
+    # Reorder arguments
+    parser.add_argument('-e', '--reorder', action='store_true', help='align the atoms of molecules (default: Hungarian)')
+    parser.add_argument('--reorder-method', action='store', default="hungarian", metavar="METHOD", help='select which reorder method to use; hungarian (default), brute, distance')
+    parser.add_argument('--use-reflections', action='store_true', help='scan through reflections in planes (eg Y transformed to -Y -> X, -Y, Z) and axis changes, (eg X and Z coords exchanged -> Z, Y, X). This will affect stereo-chemistry.')
+    parser.add_argument('--use-reflections-keep-stereo', action='store_true', help='scan through reflections in planes (eg Y transformed to -Y -> X, -Y, Z) and axis changes, (eg X and Z coords exchanged -> Z, Y, X). Stereo-chemistry will be kept.')
+
+    # Filter
+    index_group = parser.add_mutually_exclusive_group()
+    index_group.add_argument('-nh', '--no-hydrogen', action='store_true', help='ignore hydrogens when calculating RMSD')
+    index_group.add_argument('--remove-idx', nargs='+', type=int, help='index list of atoms NOT to consider', metavar='IDX')
+    index_group.add_argument('--add-idx', nargs='+', type=int, help='index list of atoms to consider', metavar='IDX')
+
+    # format and print
+    parser.add_argument('--format', action='store', help='format of input files. valid format are xyz and pdb', metavar='FMT')
+    parser.add_argument('-p', '--output', '--print', action='store_true', help='print out structure B, centered and rotated unto structure A\'s coordinates in XYZ format')
+
+    if len(sys.argv) == 1:
+        parser.print_help()
+        sys.exit(1)
+
+    args = parser.parse_args()
+
+    # As default, load the extension as format
+    if args.format is None:
+        args.format = args.structure_a.split('.')[-1]
+
+    p_all_atoms, p_all = get_coordinates(args.structure_a, args.format)
+    q_all_atoms, q_all = get_coordinates(args.structure_b, args.format)
+
+    p_size = p_all.shape[0]
+    q_size = q_all.shape[0]
+
+    if not p_size == q_size:
+        print("error: Structures not same size")
+        quit()
+
+    if np.count_nonzero(p_all_atoms != q_all_atoms) and not args.reorder:
+        msg = """
+error: Atoms are not in the same order.
+
+Use --reorder to align the atoms (can be expensive for large structures).
+
+Please see --help or documentation for more information or
+https://github.com/charnley/rmsd for further examples.
+"""
+        print(msg)
+        exit()
+
+
+    # Set local view
+    p_view = None
+    q_view = None
+
+
+    if args.no_hydrogen:
+        p_view = np.where(p_all_atoms != 'H')
+        q_view = np.where(q_all_atoms != 'H')
+
+    elif args.remove_idx:
+        index = range(p_size)
+        index = set(index) - set(args.remove_idx)
+        index = list(index)
+        p_view = index
+        q_view = index
+
+    elif args.add_idx:
+        p_view = args.add_idx
+        q_view = args.add_idx
+
+
+    # Set local view
+    if p_view is None:
+        p_coord = copy.deepcopy(p_all)
+        q_coord = copy.deepcopy(q_all)
+        p_atoms = copy.deepcopy(p_all_atoms)
+        q_atoms = copy.deepcopy(q_all_atoms)
+
+    else:
+
+        if args.reorder and args.output:
+            print("error: Cannot reorder atoms and print structure, when excluding atoms (such as --no-hydrogen)")
+            quit()
+
+        if args.use_reflections and args.output:
+            print("error: Cannot use reflections on atoms and print, when excluding atoms (such as --no-hydrogen)")
+            quit()
+
+        p_coord = copy.deepcopy(p_all[p_view])
+        q_coord = copy.deepcopy(q_all[q_view])
+        p_atoms = copy.deepcopy(p_all_atoms[p_view])
+        q_atoms = copy.deepcopy(q_all_atoms[q_view])
+
+
+    # Create the centroid of P and Q which is the geometric center of a
+    # N-dimensional region and translate P and Q onto that center.
+    # http://en.wikipedia.org/wiki/Centroid
+    p_cent = centroid(p_coord)
+    q_cent = centroid(q_coord)
+    p_coord -= p_cent
+    q_coord -= q_cent
+
+
+    # set rotation method
+    if args.rotation.lower() == "kabsch":
+        rotation_method = kabsch_rmsd
+
+    elif args.rotation.lower() == "quaternion":
+        rotation_method = quaternion_rmsd
+
+    elif args.rotation.lower() == "none":
+        rotation_method = None
+
+    else:
+        print("error: Unknown rotation method:", args.rotation)
+        quit()
+
+
+    # set reorder method
+    if not args.reorder:
+        reorder_method = None
+
+    if args.reorder_method == "hungarian":
+        reorder_method = reorder_hungarian
+
+    elif args.reorder_method == "brute":
+        reorder_method = reorder_brute
+
+    elif args.reorder_method == "distance":
+        reorder_method = reorder_distance
+
+    else:
+        print("error: Unknown reorder method:", args.reorder_method)
+        quit()
+
+
+    # Save the resulting RMSD
+    result_rmsd = None
+
+
+    if args.use_reflections:
+
+        result_rmsd, q_swap, q_reflection, q_review = check_reflections(
+            p_atoms,
+            q_atoms,
+            p_coord,
+            q_coord,
+            reorder_method=reorder_method,
+            rotation_method=rotation_method)
+
+    elif args.use_reflections_keep_stereo:
+
+        result_rmsd, q_swap, q_reflection, q_review = check_reflections(
+            p_atoms,
+            q_atoms,
+            p_coord,
+            q_coord,
+            reorder_method=reorder_method,
+            rotation_method=rotation_method,
+            keep_stereo=True)
+
+    elif args.reorder:
+
+        q_review = reorder_method(p_atoms, q_atoms, p_coord, q_coord)
+        q_coord = q_coord[q_review]
+        q_atoms = q_atoms[q_review]
+
+        if not all(p_atoms == q_atoms):
+            print("error: Structure not aligned")
+            quit()
+
+
+    # print result
+    if args.output:
+
+        if args.reorder:
+
+            if q_review.shape[0] != q_all.shape[0]:
+                print("error: Reorder length error. Full atom list needed for --print")
+                quit()
+
+            q_all = q_all[q_review]
+            q_all_atoms = q_all_atoms[q_review]
+
+        # Get rotation matrix
+        U = kabsch(q_coord, p_coord)
+
+        # recenter all atoms and rotate all atoms
+        q_all -= q_cent
+        q_all = np.dot(q_all, U)
+
+        # center q on p's original coordinates
+        q_all += p_cent
+
+        # done and done
+        xyz = set_coordinates(q_all_atoms, q_all, title="{} - modified".format(args.structure_b))
+        print(xyz)
+
+    else:
+
+        if result_rmsd:
+            pass
+
+        elif rotation_method is None:
+            result_rmsd = rmsd(p_coord, q_coord)
+
+        else:
+            result_rmsd = rotation_method(p_coord, q_coord)
+
+        print("{0}".format(result_rmsd))
+
+
+    return
+
+if __name__ == "__main__":
+    main()
diff --git a/modules/clustering.py b/modules/clustering.py
new file mode 100755
index 0000000..b0062fd
--- /dev/null
+++ b/modules/clustering.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+import sklearn
+from sklearn import cluster
+from sklearn.cluster import AffinityPropagation
+
+matrice = sys.argv[1]
+
+X = np.loadtxt(matrice) # lecture de la matrice à partir d'un fichier
+
+if len(X.shape)==0:
+    print 'Number of clusters: 1'
+    print 'Cluster centers:'
+    print '0'
+    print 'Cluster groups:'
+    print '0'
+
+else:
+    af = AffinityPropagation(affinity='precomputed').fit(X)
+
+    cluster_centers_indices = af.cluster_centers_indices_
+    labels = af.labels_
+    n_clusters_ = len(cluster_centers_indices)
+
+
+    print('Number of clusters: %d' % n_clusters_)
+    print 'Cluster centers:' 
+    print cluster_centers_indices
+
+    def ClusterIndicesNumpy(clustNum, labels_array): #numpy 
+            return np.where(labels_array == clustNum)[0]
+
+    print 'Cluster groups:'
+
+    for i in range(0,n_clusters_) :
+        a = ClusterIndicesNumpy(i, labels)
+        print a
+
+
diff --git a/modules/fe_change.sh b/modules/fe_change.sh
new file mode 100755
index 0000000..4d4b1bb
--- /dev/null
+++ b/modules/fe_change.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+molec=$1
+surf=$2
+
+clean_surf=$Surface_path/$surf.xyz
+
+surf_num_atoms=$(head -n1 $clean_surf)
+for candidate in $MolOnSurf_results_path/$molec/*; do
+	if [ -d $candidate ]; then
+		if [ ! -f $candidate/coord.xyz.old ]; then
+			cp $candidate/coord.xyz $candidate/coord.xyz.old
+		fi
+ 		cp $clean_surf $candidate/coord.xyz
+		printf "\n" >> $candidate/coord.xyz
+		tail -n +$((surf_num_atoms+3)) $candidate/coord.xyz.old >> $candidate/coord.xyz
+		molec_num_atoms=$(tail -n +$((surf_num_atoms+3)) $candidate/coord.xyz.old | wc -l)
+		sed -i "s/$surf_num_atoms/$(($surf_num_atoms+$molec_num_atoms))/g" $candidate/coord.xyz
+	fi
+done
+
diff --git a/modules/generation_conformeres.py b/modules/generation_conformeres.py
new file mode 100755
index 0000000..d76a7d8
--- /dev/null
+++ b/modules/generation_conformeres.py
@@ -0,0 +1,121 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+
+##############################################################################
+#                                                                            #
+#  Recherche de n conformeres a partir d'un fichier .mol, optimisation UFF   # 
+#                                                                            #
+#  fonctionne avec 2 arguments :                                             #
+#  - nom de la molecule (aussi nom du fichier d'entree sans .mol)            #
+#  - nombre de conformeres                                                   #
+#                                                                            #
+##############################################################################
+
+###########################
+#Generation des conformeres
+###########################
+
+
+#importation des modules nécessaires
+import sys
+import rdkit
+import os
+import glob
+from rdkit import Chem
+from rdkit.Chem import AllChem
+
+Molecule_results_path=os.environ.get('Molecule_results_path')
+Molecule_path=os.environ.get('Molecule_path')
+
+print 'Conformer generation'
+
+mol_input = sys.argv[1]  #nom du fichier d'entree
+n = int(sys.argv[2])  #nombre de conformeres
+
+m = Chem.MolFromMolFile('%s/%s.mol' % (Molecule_path, mol_input))  #lecture du fichier d'entree
+
+writer = Chem.SDWriter('%s/%s_confs.mol' % (Molecule_results_path, mol_input))
+
+m = Chem.AddHs(m)  #ajout des hydrogenes a la molecule
+m.SetProp("_Name","%s" % mol_input)  #ajout du nom de la molecule
+L = AllChem.EmbedMultipleConfs(m, numConfs=n)  #liste des identifiants des n differents conformeres
+
+for i in L:
+    AllChem.UFFOptimizeMolecule(m, confId=i)  #optimisation de la molecule pour un identifiant de conformere donne
+    writer.write(m, confId=i)  #ecriture de la geometrie du conformere dans le fichier de sortie
+
+writer.close()
+
+
+#####################################
+#Detection des conformeres identiques
+#####################################
+
+
+print 'Detection of identical conformers'
+
+liste_suppr = []
+
+for a in range(0,n): 
+    for b in range (a+1,n):
+        if a in liste_suppr:
+            break
+        else :
+            if a != b :
+                valeur = AllChem.GetBestRMS(m, m, prbId=b, refId=a) 
+                if valeur < 0.01: #valeur de RMSD limite en dessous de laquelle les conformeres sont consideres comme identiques
+                    liste_suppr.append(a)
+
+
+#################################
+#Separation en plusieurs fichiers 
+#################################
+
+
+#fait uniquement pour les conformeres differents
+
+if not os.path.exists('%s/%s_confs' % (Molecule_results_path, mol_input)):  #verification de l'existance du dossier dans lequel les n fichiers seront ecrits
+    os.mkdir('%s/%s_confs' % (Molecule_results_path, mol_input))  #creation du dossier si inexistant
+else :
+    liste=glob.glob('%s/%s_confs/*.mol' % (Molecule_results_path, mol_input))
+    for file in liste:
+        os.remove(file)
+
+
+i = 1
+liste_ok = []
+liste_non_ok = []
+
+nb_lines = sum(1 for line in open('%s/%s_confs.mol' % (Molecule_results_path, mol_input))) #lecture du nombre de lignes total du fichier contenant les n conformeres
+lg = nb_lines/n #calcul du nombre de lignes pour un conformere
+
+for j in range(0,n):
+    if j not in liste_suppr :
+        liste_ok.append(j)
+        l = 1 #compteur du numero de ligne
+        lg_init = 1 + lg * j  #numero de la ligne a laquelle la structure du conformere commence
+        lg_fin = lg * (j + 1) #numero de la ligne a laquelle la structure du conformere termine
+        line = open('%s/%s_confs.mol' % (Molecule_results_path, mol_input)).readline()
+        with open('%s/%s_confs.mol' % (Molecule_results_path, mol_input)) as fp:
+            #si la ligne est comprise entre la ligne a laquelle la structure du conformere commence et la ligne a laquelle elle termine alors cette ligne est copiee dans le fichier correspondant au numero du conformere
+            for line in fp :
+                if l <= lg_fin :
+                    if l >= lg_init :
+                        file = open('%s/%s_confs/%s_%s.mol' % (Molecule_results_path, mol_input, mol_input, i),"a")
+                        file.write(line)
+                        file.close()
+                    l += 1
+                else :
+                    break
+        i += 1
+    else:
+        liste_non_ok.append(j)
+
+os.remove('%s/%s_confs.mol' % (Molecule_results_path, mol_input))
+
+
+print ('number of removed conformers : %s' % (len(liste_non_ok)))
+#print liste_non_ok
+print ('number of conformers considered as different : %s' % (len(liste_ok)))
+
diff --git a/modules/launch_cp2k_molecule_seule.sh b/modules/launch_cp2k_molecule_seule.sh
new file mode 100755
index 0000000..10ec7a4
--- /dev/null
+++ b/modules/launch_cp2k_molecule_seule.sh
@@ -0,0 +1,31 @@
+#!/bin/bash
+
+molecule=$1
+
+n="$(find ${Molecule_results_path}/${molecule}_confs/xyz_files/${molecule}*.xyz | wc -l)"
+
+if [ ! -d "${Molecule_results_path}/${molecule}" ]
+then mkdir ${Molecule_results_path}/${molecule}
+else 
+	print 'problem there is already a directory for this molecule'
+	exit
+fi
+
+mkdir ${Molecule_results_path}/${molecule}/${molecule}_1
+cd ${Molecule_results_path}/${molecule}/${molecule}_1
+sed "/PROJECT_NAME/c\PROJECT_NAME ${molecule}" ${CP2K_input_molecule} > ${Molecule_results_path}/${molecule}/${molecule}_1/${molecule}.inp
+sed "s/XXXX/${molecule}/g" ${CP2K_sub} > ${Molecule_results_path}/${molecule}/${molecule}_1/cp2k_gamma.j
+cp ${Molecule_results_path}/${molecule}_confs/xyz_files/${molecule}_1.xyz ${Molecule_results_path}/${molecule}/${molecule}_1/coord.xyz
+qsub cp2k_gamma.j
+
+for ((i=2; i<=n; i++)) ; do
+	mkdir ${Molecule_results_path}/${molecule}/${molecule}_$i
+	cd ${Molecule_results_path}/${molecule}/${molecule}_$i
+	ln ${Molecule_results_path}/${molecule}/${molecule}_1/cp2k_gamma.j .
+	ln ${Molecule_results_path}/${molecule}/${molecule}_1/${molecule}.inp .
+	cp ${Molecule_results_path}/${molecule}_confs/xyz_files/${molecule}_${i}.xyz coord.xyz
+	qsub cp2k_gamma.j
+done
+
+
+
diff --git a/modules/launch_script+diss.sh b/modules/launch_script+diss.sh
new file mode 100755
index 0000000..f412c35
--- /dev/null
+++ b/modules/launch_script+diss.sh
@@ -0,0 +1,48 @@
+#!/bin/bash
+
+nom_de_la_molecule=$1
+
+mol_dir=${MolOnSurf_results_path}/${nom_de_la_molecule}
+
+list_errors=()
+
+if [ -e ${mol_dir}/errors ]
+then
+	nb_lines="$(cat ${mol_dir}/errors | wc -l)"
+	for ((line=1; line<=$nb_lines; line++ )) ; do
+		err_dir="$(awk -v var="$line" 'NR == var' ${mol_dir}/errors | awk '{print $2}')"
+		list_errors+=(${err_dir})
+	done
+fi
+
+lg_liste=${#list_errors[@]}
+
+sed "s/XXXX/surf_${nom_de_la_molecule}/g" ${CP2K_sub} > ${mol_dir}/${nom_de_la_molecule}_1/cp2k_gamma.j
+sed "s/molecule/${nom_de_la_molecule}/g" ${CP2K_input_MolOnSurf} > ${mol_dir}/${nom_de_la_molecule}_1/surf_${nom_de_la_molecule}.inp
+
+for dir in ${mol_dir}/${nom_de_la_molecule}* ; do
+	i_error=0
+	if [ -e ${mol_dir}/errors ]
+	then 
+		for ((j=0; j<$lg_liste; j++)) ; do
+			if [[ "$dir" == *"${list_errors[$j]}" ]]
+			then 
+				i_error=1 
+				echo "error"
+			fi
+		done
+	fi
+	if [ "$i_error" -eq "0" ]
+	then 
+		cd $dir
+		if [ "${dir}" != "${mol_dir}/${nom_de_la_molecule}_1"  ]
+		then 
+			ln ${mol_dir}/${nom_de_la_molecule}_1/cp2k_gamma.j .
+			ln ${mol_dir}/${nom_de_la_molecule}_1/surf_${nom_de_la_molecule}.inp .
+		fi
+		qsub cp2k_gamma.j
+	fi
+done
+
+
+	
diff --git a/modules/mol_to_xyz.sh b/modules/mol_to_xyz.sh
new file mode 100755
index 0000000..98a1273
--- /dev/null
+++ b/modules/mol_to_xyz.sh
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+molecule=$1
+
+n="$(find ${Molecule_results_path}/${molecule}_confs/${molecule}* | wc -l)"
+
+ligne="$(awk 'NR==4' ${Molecule_results_path}/${molecule}_confs/${molecule}_1.mol)"
+nb_atomes="$(echo $ligne | awk '{print $1}')"
+nb_ligne_fin=$nb_atomes
+let nb_ligne_fin+=4
+ 
+if [ ! -d "${Molecule_results_path}/${molecule}_confs/xyz_files" ]
+	then mkdir ${Molecule_results_path}/${molecule}_confs/xyz_files
+	else rm ${Molecule_results_path}/${molecule}_confs/xyz_files/*
+	fi
+
+for ((j=1; j<=n; j++)) ; do
+	mol_file=${Molecule_results_path}/${molecule}_confs/${molecule}_$j.mol
+	xyz_file=${Molecule_results_path}/${molecule}_confs/xyz_files/${molecule}_$j.xyz
+
+	for ((i=5; i<=nb_ligne_fin; i++)) ; do
+		awk 'NR==n' n=$i $mol_file >> result
+	done
+
+	awk '{print $4, " ", "/", $1, "/" , "/" , $2, "/" , "/" , $3, "/" }' result >> result3
+	
+	echo $nb_atomes >> $xyz_file
+	echo ' ' >> $xyz_file
+
+	column -s / -t result3 >> $xyz_file
+
+	rm result
+	rm result3
+done
+
diff --git a/modules/recenter_only.py b/modules/recenter_only.py
new file mode 100755
index 0000000..8cdff06
--- /dev/null
+++ b/modules/recenter_only.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+import os
+import ase
+from ase import io
+from ase.io import read
+from ase.io import write
+from ase import Atom
+from ase import build
+from ase import geometry
+from ase.geometry import wrap_positions
+import ase_Fe1_Fe2
+
+
+molsurf_dir = sys.argv[1]
+central_atom_number = int(sys.argv[2]) # atome sur lequel la surface+molecule va etre recentrée.
+molsurf_file = molsurf_dir + '/last_geo.xyz'
+recentered_file = molsurf_dir + '/recentered_only.xyz'
+
+molsurf = read(molsurf_file)
+vect_translation = - molsurf[central_atom_number].position
+vect_translation[2] = 0
+molsurf.translate(vect_translation)
+
+write(recentered_file , molsurf)
+
+
diff --git a/modules/recenter_periodicity.py b/modules/recenter_periodicity.py
new file mode 100755
index 0000000..7abc052
--- /dev/null
+++ b/modules/recenter_periodicity.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+import os
+import ase
+from ase import io
+from ase.io import read
+from ase.io import write
+from ase import Atom
+from ase import build
+from ase import geometry
+from ase.geometry import wrap_positions
+import ase_Fe1_Fe2
+
+a1=float(os.environ.get('a1'))
+a2=float(os.environ.get('a2'))
+a3=float(os.environ.get('a3'))
+b1=float(os.environ.get('b1'))
+b2=float(os.environ.get('b2'))
+b3=float(os.environ.get('b3'))
+c1=float(os.environ.get('c1'))
+c2=float(os.environ.get('c2'))
+c3=float(os.environ.get('c3'))
+
+molsurf_dir = sys.argv[1]
+central_atom_number = int(sys.argv[2]) # atome sur lequel la surface+molecule va etre recentrée.
+
+molsurf_file = molsurf_dir + '/last_geo.xyz'
+recentered_file = molsurf_dir + '/recentered_periodicity.xyz'
+
+molsurf = read(molsurf_file)
+vect_translation = - molsurf[central_atom_number].position
+vect_translation[2] = 0
+molsurf.translate(vect_translation)
+
+molsurf.set_pbc([1,1,0])
+molsurf.set_cell([[a1,a2,a3],[b1,b2,b3],[c1,c2,c3]])
+
+molsurf.wrap()
+
+write(recentered_file , molsurf)
+
+
+
+
+
+
+
+
+
+
diff --git a/modules/script_add_adsorbate+diss.sh b/modules/script_add_adsorbate+diss.sh
new file mode 100755
index 0000000..b89d20a
--- /dev/null
+++ b/modules/script_add_adsorbate+diss.sh
@@ -0,0 +1,83 @@
+#!/bin/bash
+
+# script surface.xyz at.surf molecule.xyz at.mol distance at.surf2 at.mol2 alpha theta output(.xyz)
+
+distance_mol_surf=2
+
+nom_de_la_molecule=$1
+surface_d_adsorption=$2
+list_atom_1_mol=($3)
+list_atom_2_mol=($4)
+list_atom_surf=($5)
+
+nombre_atom_mol=${#list_atom_1_mol[@]}
+
+angles_phi=(0 0 0 270 270 270 180 180 180 90 90 90)
+angles_theta=(0 300 300 0 300 300 0 300 300 0 300 300)
+angles_psi=(0 0 315 0 0 315 0 0 315 0 0 315)
+
+nombre_atom_surf=${#list_atom_surf[@]}
+nombre_angles=${#angles_phi[@]}
+
+if [ ! -d ${MolOnSurf_results_path}/${nom_de_la_molecule} ]
+then mkdir ${MolOnSurf_results_path}/${nom_de_la_molecule}
+else
+	echo "There is already a directory for this molecule: remove this directory?"
+	read reponse
+	if [[ "$reponse" == "yes" ]]
+	then rm -r ${MolOnSurf_results_path}/${nom_de_la_molecule}/*
+	fi 
+fi
+
+
+##################
+# Détermination du numéro correspondant à la molécule ayant la structure la plus stable
+##################
+
+
+n="$(find ${Molecule_results_path}/${nom_de_la_molecule}/${nom_de_la_molecule}*/coord.xyz | wc -l)"
+min_energy_molecule=0
+
+for ((a=1; a<=n; a++)) ; do
+	energy_molecule="$(awk '/ENERGY/{E=$NF} END{print E}' ${Molecule_results_path}/${nom_de_la_molecule}/${nom_de_la_molecule}_$a/*${nom_de_la_molecule}.out)"
+	if [[ "$energy_molecule" > "$min_energy_molecule" ]]
+	then 
+		min_energy_molecule=$energy_molecule
+		numero_mol_stable=$a
+	fi
+done
+
+
+##################
+# Génération des structures de molécules adsorbées sur la surface
+##################
+
+b=1
+
+mkdir ${MolOnSurf_results_path}/${nom_de_la_molecule}/coords
+
+molecule_a_adsorber=${nom_de_la_molecule}_${numero_mol_stable}
+for ((i=0; i<${nombre_atom_surf}; i++)) ; do
+	atom_surf=${list_atom_surf[$i]}
+	for ((j=0; j<${nombre_atom_mol}; j++)) ; do
+		atom_1_mol=${list_atom_1_mol[$j]}
+		atom_2_mol=${list_atom_2_mol[$j]}
+		echo " "
+		echo "Adsorption of atom $atom_1_mol on atom $atom_surf" 
+		for ((k=0; k<${nombre_angles}; k++)) ; do
+			angle_phi=${angles_phi[$k]}
+			angle_theta=${angles_theta[$k]}
+			angle_psi=${angles_psi[$k]}
+			${DockOnSurf_path}/modules/add_adsorbate_euler+dissociation.py ${Surface_path}/${surface_d_adsorption}.xyz ${atom_surf} ${Molecule_results_path}/${nom_de_la_molecule}/${molecule_a_adsorber}/coord.xyz ${atom_1_mol} ${distance_mol_surf} ${atom_2_mol} ${angle_phi} ${angle_theta} ${angle_psi} ${MolOnSurf_results_path}/${nom_de_la_molecule}/coords/${nom_de_la_molecule}_$b
+			b=$((b+1))
+		done
+	done
+done
+
+for file_molsurf in ${MolOnSurf_results_path}/${nom_de_la_molecule}/coords/*; do
+	name="$(echo $file_molsurf | awk -F/ '{print $NF}' | awk -F. '{print $(NF-1)}')"
+	mkdir ${MolOnSurf_results_path}/${nom_de_la_molecule}/$name
+	mv $file_molsurf ${MolOnSurf_results_path}/${nom_de_la_molecule}/$name/coord.xyz
+done
+rm -r ${MolOnSurf_results_path}/${nom_de_la_molecule}/coords
+
diff --git a/modules/script_grandes_molecules+diss.sh b/modules/script_grandes_molecules+diss.sh
new file mode 100755
index 0000000..9f8d69e
--- /dev/null
+++ b/modules/script_grandes_molecules+diss.sh
@@ -0,0 +1,363 @@
+#!/bin/bash
+
+##########################
+#### Arguments input  ####
+##########################
+
+nom_de_la_molecule=$1
+atom_nb=$2
+nb_surf=$3
+cutoff=$4
+
+
+#############################################################
+#### Creation of the working diretory and initialisation ####
+#############################################################
+
+
+mol_dir=${MolOnSurf_results_path}/${nom_de_la_molecule}
+
+list_errors=()
+
+mkdir ${mol_dir}/analyse
+
+
+###############################################################
+#### Extraction of structures needed from file *-pos-1.xyz ####
+###############################################################
+
+
+num=0
+
+for dir in ${mol_dir}/* ; do
+	if [ ! -e  $dir/*-pos-1.xyz ]
+	then list_errors+=($dir)
+	else
+		name="$(echo $dir | awk -F/ '{print $NF}')"
+		mkdir ${mol_dir}/analyse/${name}_$num
+
+
+		#######################################
+		#### Get last geometry of mol+surf #### 
+		#######################################
+
+
+		n_line="$(grep -n "i =" ${dir}/*-pos-1.xyz | awk 'END{print}' | awk 'NF=1{print}' | awk 'sub(".$", "")')"
+		n_line_1="$(expr $n_line - 1)"
+		awk -v ligne="$n_line_1" 'NR>=ligne{print}' ${dir}/*-pos-1.xyz > ${mol_dir}/analyse/${name}_$num/last_geo.xyz
+
+
+		#######################################################
+		#### Recenter the mol+surf on one atom without pbc ####
+		#######################################################
+
+		
+		${DockOnSurf_path}/modules/recenter_only.py ${mol_dir}/analyse/${name}_$num ${atom_nb}
+		recentered_temp_file=${mol_dir}/analyse/${name}_$num/recentered_only.xyz
+
+
+		#####################################################
+		#### Get the molecule alone centered on one atom ####
+		#####################################################
+
+		nb_surf2="$(expr ${nb_surf} + 2 )"
+		nbligne="$(awk 'END{print NR}' ${recentered_temp_file})"
+		for ((a=3; a<=nbligne; a++)) ; do
+			z="$(awk '{print NR,$0}' ${recentered_temp_file}| awk 'NR == nligne {print}' nligne="$a" | awk '{print $1}')"
+			if (( $(echo "$z > ${nb_surf2}" | bc -l) ))
+			then awk 'NR == nligne {print}' nligne="$a" ${recentered_temp_file} >> ${mol_dir}/analyse/${name}_$num/temp_file.xyz
+			fi
+		done
+
+       		temp_file=${mol_dir}/analyse/${name}_$num/temp_file.xyz
+		awk 'END{print NR}' ${temp_file} > ${mol_dir}/analyse/${name}_$num/molecule_seule.xyz
+       		mol_seule_file=${mol_dir}/analyse/${name}_$num/molecule_seule.xyz
+		echo ' ' >> ${mol_seule_file}
+		cat ${temp_file} >> ${mol_seule_file}
+		rm ${temp_file}
+
+
+		##############################################################
+		### Apply pbc conditions to surface only in mol+surf file ####
+		##############################################################
+
+
+		${DockOnSurf_path}/modules/recenter_periodicity.py ${mol_dir}/analyse/${name}_$num ${atom_nb}
+		recentered_periodicity_file=${mol_dir}/analyse/${name}_$num/recentered_periodicity.xyz
+
+		for ((a=3; a<=nbligne; a++)) ; do
+			z="$(awk '{print NR,$0}' ${recentered_temp_file}| awk 'NR == nligne {print}' nligne="$a" | awk '{print $1}')"
+			if (( $(echo "$z <= ${nb_surf2}" | bc -l) ))
+			then awk 'NR == nligne {print}' nligne="$a" ${recentered_periodicity_file} >> ${mol_dir}/analyse/${name}_$num/temp_file.xyz
+			fi
+		done
+
+		awk 'NR <= 2{print}' ${recentered_periodicity_file} > ${mol_dir}/analyse/${name}_$num/recentered_pbc.xyz
+		cat ${mol_dir}/analyse/${name}_$num/temp_file.xyz >> ${mol_dir}/analyse/${name}_$num/recentered_pbc.xyz
+		awk 'NR > 2{print}' ${mol_seule_file} >> ${mol_dir}/analyse/${name}_$num/recentered_pbc.xyz
+
+		rm ${mol_dir}/analyse/${name}_$num/temp_file.xyz
+		rm ${recentered_periodicity_file}
+		rm ${recentered_temp_file}
+
+		
+		num=$((num+1))
+
+	fi
+done
+
+echo "All structures files have been created"
+
+
+#################################################################################
+#### RMSD calculation of molecules without surface to create the RMSD matrix ####
+#################################################################################
+
+
+matrix_file=${mol_dir}/analyse/matrice_RMSD_mol_seule.txt
+
+touch ${matrix_file}
+
+for ((i=0; i<$num; i++)) ; do
+	echo  >> ${matrix_file}
+done
+
+for ((i=1; i<=$num; i++)) ; do
+	for ((j=$i; j<=$num; j++)) ; do
+		if (( $(echo "$i == $j" | bc -l) ))
+		then 
+			add_to_line="$(awk -v ii="$i" 'NR==ii {print}' ${matrix_file})"
+			add_to_line+=' 0'
+			sed -i "${i}c $add_to_line" ${matrix_file}
+		else 
+			num_i="$(expr $i - 1)"
+			num_j="$(expr $j - 1)"
+			RMSD_mol_seule_i_j="$(${DockOnSurf_path}/modules/calculate_rmsd --reorder ${mol_dir}/analyse/${nom_de_la_molecule}*_${num_i}/molecule_seule.xyz ${mol_dir}/analyse/${nom_de_la_molecule}*_${num_j}/molecule_seule.xyz)"
+			add_to_i_line="$(awk -v ii="$i" 'NR==ii {print}' ${matrix_file})"
+			add_to_i_line+=' -'
+			add_to_i_line+=$RMSD_mol_seule_i_j
+			sed -i "${i}c $add_to_i_line" ${matrix_file}
+			add_to_j_line="$(awk -v jj="$j" 'NR==jj {print}' ${matrix_file})"
+			add_to_j_line+=' -'
+			add_to_j_line+=$RMSD_mol_seule_i_j
+			sed -i "${j}c $add_to_j_line" ${matrix_file}
+		fi
+	done
+done
+
+echo "RMSD matrix for molecules without surface has been calculated"
+
+
+###############################################################
+#### Clustering of the molecule structures without surface #### 
+###############################################################
+
+
+cluster_molecule_file=${mol_dir}/analyse/cluster_molecule.txt
+
+${DockOnSurf_path}/modules/clustering.py ${matrix_file} > ${cluster_molecule_file}
+
+echo "Clustering of molecules has been done"
+
+
+###################################################################################
+#### Separation of the clusters found for the molecule alone and re-clustering ####
+###################################################################################
+
+
+nb_clusters="$(awk 'END{print NR}' ${cluster_molecule_file})"
+
+sed -i 's/\[/ /g' ${cluster_molecule_file}
+sed -i 's/\]/ /g' ${cluster_molecule_file}
+
+mkdir ${mol_dir}/analyse/cluster_centers_energy
+
+s=1
+
+line_cluster_group="$(awk '{print NR,$0}' ${cluster_molecule_file} | grep "Cluster\ groups:" | awk '{print $1}')"
+line_first_cluster="$(expr ${line_cluster_group} + 1 )"
+
+for ((k=${line_first_cluster}; k<=${nb_clusters}; k++)) ; do
+
+
+	################################
+	#### Separation of clusters ####
+	################################
+
+
+	awk -v ligne="$k" 'NR==ligne{print}' ${cluster_molecule_file} > ${mol_dir}/analyse/array
+	sed -i 's/^ *//' ${mol_dir}/analyse/array
+	sed -i 's/  / /g' ${mol_dir}/analyse/array
+	sed -i 's/  / /g' ${mol_dir}/analyse/array
+	sed -i 's/  / /g' ${mol_dir}/analyse/array
+	sed -i 's/  / /g' ${mol_dir}/analyse/array
+	readarray -d " " cluster < ${mol_dir}/analyse/array
+	
+	mkdir ${mol_dir}/analyse/cluster_molecule_$s
+	cluster_dir=${mol_dir}/analyse/cluster_molecule_$s
+	
+	for numero in ${cluster[*]} ; do
+		mv ${mol_dir}/analyse/${nom_de_la_molecule}*_${numero} ${cluster_dir}/.
+	done
+	a=0
+	for dir2 in ${cluster_dir}/${nom_de_la_molecule}* ; do
+		mv $dir2 ${dir2}_$a
+		a=$((a+1))
+	done
+	
+	rm ${mol_dir}/analyse/array
+
+
+	######################################################################
+	#### Calculation of RMSD for mol+surf and creation of RMSD matrix ####
+	######################################################################
+
+	
+ 	cluster_matrix_file=${cluster_dir}/matrice_RMSD_cluster.txt
+ 	touch ${cluster_matrix_file}
+
+	len=${#cluster[@]}
+	len_cluster="$(expr $len - 1)"
+
+ 	for ((i=1; i<=${len_cluster}; i++)) ; do
+ 		echo >> ${cluster_matrix_file}
+ 	done
+
+ 	for ((i=1; i<=${len_cluster}; i++)) ; do
+ 		for ((j=$i; j<=${len_cluster}; j++)) ; do
+ 			if (( $(echo "$i == $j" | bc -l) ))
+ 			then
+ 				add_to_line="$(awk -v ii="$i" 'NR==ii {print}' ${cluster_matrix_file})"
+ 				add_to_line+=' 0'
+ 				sed -i "${i}c $add_to_line" ${cluster_matrix_file}
+ 			else
+				num_i="$(expr $i - 1)"
+				num_j="$(expr $j - 1)"
+ 				RMSD_mol_surf_i_j="$(${DockOnSurf_path}/modules/calculate_rmsd --reorder ${cluster_dir}/${nom_de_la_molecule}*_${num_i}/recentered_pbc.xyz ${cluster_dir}/${nom_de_la_molecule}*_${num_j}/recentered_pbc.xyz)"
+ 				add_to_i_line="$(awk -v ii="$i" 'NR==ii {print}' ${cluster_matrix_file})"
+ 				add_to_i_line+=' -'
+ 				add_to_i_line+=$RMSD_mol_surf_i_j
+ 				sed -i "${i}c $add_to_i_line" ${cluster_matrix_file}
+ 				add_to_j_line="$(awk -v jj="$j" 'NR==jj {print}' ${cluster_matrix_file})"
+ 				add_to_j_line+=' -'
+ 				add_to_j_line+=$RMSD_mol_surf_i_j
+ 				sed -i "${j}c $add_to_j_line" ${cluster_matrix_file}
+ 			fi
+ 		done
+ 	done
+
+
+	#######################
+	#### Re-clustering ####
+	#######################
+
+	
+	cluster_molsurf_file=${cluster_dir}/cluster_molsurf.txt
+	${DockOnSurf_path}/modules/clustering.py ${cluster_matrix_file} > ${cluster_molsurf_file}
+
+	
+	#################################################
+	#### Get the cluster centers (min of energy) ####
+	#################################################
+
+	
+	nb_cluster_molsurf="$(awk 'END{print NR}' ${cluster_molsurf_file})"
+	line_cluster_group2="$(awk '{print NR,$0}' ${cluster_molsurf_file} | grep "Cluster\ groups:" | awk '{print $1}')"
+	line_first_cluster2="$(expr ${line_cluster_group2} + 1 )"
+
+	for ((a=${line_first_cluster2}; a<=${nb_cluster_molsurf}; a++)); do
+		awk -v ligne="$a" 'NR==ligne{print}' ${cluster_molsurf_file} > ${cluster_dir}/array
+		sed -i 's/\[/ /g' ${cluster_dir}/array
+		sed -i 's/\]/ /g' ${cluster_dir}/array
+		sed -i 's/^ *//' ${cluster_dir}/array
+		sed -i 's/  / /g' ${cluster_dir}/array
+		readarray -d " " cluster_struct < ${cluster_dir}/array
+	
+		E_min=0
+		numero_mol_E_min=t
+
+		for struct_nb in ${cluster_struct[*]} ; do
+			energie="$(awk 'NR==2{E=$NF} END{print E}' ${cluster_dir}/*_${struct_nb}/last_geo.xyz)"
+			if (( $(echo $energie '<' $E_min | bc -l) ))
+			then
+				E_min=$energie
+				numero_mol_E_min=$struct_nb
+			fi
+		done
+
+		cp -r ${cluster_dir}/${nom_de_la_molecule}*_${numero_mol_E_min} ${mol_dir}/analyse/cluster_centers_energy/.
+	done
+
+	s=$((s+1))
+done
+
+for dir in ${mol_dir}/analyse/cluster_centers_energy/${nom_de_la_molecule}* ; do
+	nb_champs="$(echo $dir | awk -F  "_" '{print NF}')"
+	nb_limit="$(expr $nb_champs - 2)"
+	name="$(echo $dir | cut -d '_' -f-$nb_limit)"
+	mv $dir $name
+done
+
+echo "Clustering of structures mol+surf has been done"
+
+
+###################################################
+#### Determination of calculations to relaunch ####
+###################################################
+
+Energy_min=0
+
+for structure in ${mol_dir}/analyse/cluster_centers_energy/*; do
+	Energy="$(awk 'NR==2{E=$NF} END{print E}' $structure/last_geo.xyz)"
+	if (( $(echo $Energy '<' $Energy_min | bc -l) ))
+	then Energy_min=$Energy
+	fi
+done
+
+Limit_energy="$(echo $Energy_min + ${cutoff}/27.2 | bc -l)" #limit of energy to relaunch the calculation 0.25eV more than min
+
+mkdir ${mol_dir}/relaunched_calculations
+
+for structure in ${mol_dir}/analyse/cluster_centers_energy/*; do
+	Energy="$(awk 'NR==2{E=$NF} END{print E}' $structure/last_geo.xyz)"
+	if (( $(echo $Energy '<' $Limit_energy | bc -l) ))
+	then cp -r $structure ${mol_dir}/relaunched_calculations/.
+	fi
+done
+
+
+####################################
+#### Launching the calculations ####
+####################################
+
+n="$(find ${mol_dir}/relaunched_calculations/${nom_de_la_molecule}* -type d | wc -l)"
+
+echo $n "clusters centers with energy lower than Emin+$cutoff eV have been found"
+echo "Relaunch the calculations ?"
+read answer
+
+if [[ $answer = "yes" ]]
+then 
+	cp ${mol_dir}/${nom_de_la_molecule}_1/cp2k_gamma.j ${mol_dir}/relaunched_calculations/.
+
+	for dir in ${mol_dir}/relaunched_calculations/${nom_de_la_molecule}*; do
+                # takes the restart file of the previous calculation as input file for re-optimization
+		conf=$(echo $dir | rev | cut -d/ -f1 | rev )
+                cp ${mol_dir}/$conf/surf_${nom_de_la_molecule}-1.restart ${dir}/surf_${nom_de_la_molecule}.inp
+                topo_line=$(grep -n \&TOPOLOGY ${dir}/surf_${nom_de_la_molecule}.inp | cut -d: -f1)
+                sed -i "$topo_line,$(($topo_line+5))d" ${dir}/surf_${nom_de_la_molecule}.inp
+		# changes the optimization criterions for more precise re-optimization
+		sed -i 's/MAX_FORCE/#MAX_FORCE/g' ${dir}/surf_${nom_de_la_molecule}.inp
+		sed -i 's/MAX_DR/#MAX_DR/g' ${dir}/surf_${nom_de_la_molecule}.inp
+		sed -i 's/RMS_DR/#RMS_DR/g' ${dir}/surf_${nom_de_la_molecule}.inp
+		sed -i 's/RMS_FORCE/#RMS_FORCE/g' ${dir}/surf_${nom_de_la_molecule}.inp
+		sed -i 's/EPS_SCF/EPS_SCF\ 1E-07\ #/g' ${dir}/surf_${nom_de_la_molecule}.inp
+		# launches the calculation
+		cd $dir/
+		qsub ${mol_dir}/relaunched_calculations/cp2k_gamma.j
+		cd -
+	done
+fi
+
+
+
