First make 1 cluster:

~/md/structure/makeFCC 3.62 12 | grep ^Cu | awk '{ r=sqrt($2^2+$3^2+$4^2); if (r<10) print $2,$3,$4 }' > clus_10A

m clus_10A 
-8.145 -4.525 -2.715
-8.145 -4.525 0.905
-8.145 -2.715 -4.525
-8.145 -0.905 -2.715
-8.145 -2.715 -0.905
-8.145 -0.905 0.905

  Then rotate twice to get two different clusters (script given below):

rotate_euler 34 97 22 < clus_10A > clus_10A_rot1
rotate_euler 24 17 62 < clus_10A > clus_10A_rot2

  The form joint xyz file for mdmorse:

wc clus_10A_rot1
    336    1008    9336 clus_10A_rot1
wc clus_10A_rot2
    336    1008    9336 clus_10A_rot2

  Check what is actual minimum and maximum z coordinate:

min < clus_10A_rot1                             
-9.08003 -9.16872 -9.32828 
max < clus_10A_rot2 
9.33743 9.70212 9.32001 

  Using this information placed two clusters so they are separated by 4 Å in z:
This means old minimum or maximum of about 9.3 should be displaced
to +2 and -2 Å, i.e. by 11.3 Å:

cat clus_10A_rot1 clus_10A_rot2 | awk 'BEGIN { print 2*336; print } { if (NR<=336) print "Cu",$1,$2,$3+11.3,1; else print "Cu",$1,$2,$3-11.3,1; }' > atoms.in

  Checked by rasmol this is OK.

  Then started simulation.

  Used no P control, btctau=300, open boundaries in all dimensions.

 mdmorse |& tee out.600K_1

   Used cutoff=7.0 since this is pretty close to accurate (compare
previous simulations on course)

dpc msleep 100 m 4 d 20 cube 20 sd 800 800 xyz erase sort 2 4 3 5 atoms.out

  Answer: the two clusters merge and the lattice planes join.
Actually this happens many times, the cluster transforms shape
spontaneously!

  Then restarted with extremely slow cooling:

tail -674 < atoms.out | awk 'NF<3 { print; } NF>2 { print $1,$2,$3,$4,1 }' > atoms.in

  Set btctau to 30000, then
 
mdmorse | & tee out.600K_1_cool


-----------------------------------------------------------------------------------------------------

  Script rotate_euler:


gawk -v psi=$1 -v phi=$2 -v theta=$3 'BEGIN { 

	if (psi==-9999 && phi=-9999 && theta=-9999) {
		srand();
		psi=rand()*180.0;
		phi=rand()*180.0;
		theta=rand()*180.0;
		printf "rotate_euler using random angle %g %g %g\n",psi,phi,theta >> "/dev/stderr";
	}

	psi=3.1415927/180*phi; 
	phi=3.1415927/180*psi; 
	theta=3.1415927/180*theta; 
} 
{ 
x=$1*(cos(phi)*cos(psi)-sin(phi)*sin(psi)*cos(theta))+$2*(-sin(phi)*cos(psi)-cos(phi)*sin(psi)*cos(theta))+$3*sin(psi)*sin(theta);

y=$1*(cos(phi)*sin(psi)+sin(phi)*cos(psi)*cos(theta))+$2*(-sin(phi)*sin(psi)+cos(phi)*cos(psi)*cos(theta))-$3*cos(psi)*sin(theta);

z=$1*sin(phi)*sin(theta)+$2*cos(phi)*sin(theta)+$3*cos(theta);

printf "%g %g %g %s %s\n",x,y,z,$4,$5,$6
}'