Chimie Théorique » scripts_chimie4psmn » DockOnSurf

#!/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/

		new_job_id=$(qsub ${mol_dir}/relaunched_calculations/cp2k_gamma.j | awk '{print $3}')

    echo "Submitted job $new_job_id"

    jobs_owned+=($new_job_id)

		cd -

	done

fi

### Attente jusqu'à ce que tous les calculs aient fini

go_on=true

while [ $go_on == true ]; do

  all_jobs=`qstat | tail -n+3 | awk '{print $1}'`

  for j1 in ${jobs_owned[@]} ; do

    for j2 in ${all_jobs[@]}; do

      if [ $j1 == $j2 ]; then

        sleep 30

        continue 3

      fi

    done

  done

  echo "loop finished"

  go_on=false

done