non-orthogonal PBC

Author: slipchenko

.idea/.gitignore

@@ -0,0 +1,77 @@
+#! /usr/bin/perl
+
+# This script generates a rectangular box of fragments for EFPMD input
+
+use 5.008;
+use strict;
+use warnings;
+
+use Math::Trig;
+
+if (scalar(@ARGV) != 6) {
+	die <<EOF;
+usage: cubegen.pl <n1:n2:...> <r1:r2:...> <space> <nx> <ny> <nz>
+
+  <n1:n2:...>  a colon separated list of fragment names
+  <r1:r2:...>  a colon separated list of fragment ratios
+      <space>  distance between the fragments
+         <nx>  number of fragments in x direction
+         <ny>  number of fragments in y direction
+         <nz>  number of fragments in z direction
+
+example: cubegen.pl c2h5oh_l:h2o_l 40:60 4.0 20 20 20 > vodka.in
+EOF
+}
+
+my @name = split ':', $ARGV[0];
+my @ratio = split ':', $ARGV[1];
+my $space = $ARGV[2];
+my $nx = $ARGV[3];
+my $ny = $ARGV[4];
+my $nz = $ARGV[5];
+
+die "ratios do not match names" unless scalar(@name) == scalar(@ratio);
+die "positive number expected" unless $space > 0.0 && $nx > 0 && $ny > 0 && $nz > 0;
+
+foreach my $i (1 .. scalar(@ratio) - 1) {
+	$ratio[$i] += $ratio[$i - 1];
+}
+
+foreach my $x (0 .. $nx - 1) {
+	foreach my $y (0 .. $ny - 1) {
+		foreach my $z (0 .. $nz - 1) {
+			print_fragment($x, $y, $z);
+		}
+	}
+}
+
+sub select_fragment {
+	my $n = int(rand($ratio[-1]));
+
+	foreach my $i (0 .. scalar(@ratio) - 1) {
+		if ($ratio[$i] > $n) {
+			return $i;
+		}
+	}
+
+	die;
+}
+
+sub print_fragment {
+	my ($x, $y, $z) = @_;
+	my $idx = select_fragment();
+	my @xyzabc = ($space * $x,
+		      $space * $y,
+		      $space * $z,
+		      2.0 * pi * rand,
+		      pi * rand,
+		      2.0 * pi * rand);
+
+	print "\n", "fragment ", $name[$idx], "\n";
+
+	foreach (@xyzabc) {
+		printf " %8.3f", $_;
+	}
+
+	print "\n";
+}

.idea/libefp_pbc.iml

@@ -0,0 +1,31 @@
+#! /usr/bin/perl
+
+# This script extracts trajectory data from EFPMD output
+
+use 5.008;
+use strict;
+use warnings;
+
+die "usage: trajectory.pl <input>\n" if scalar(@ARGV) != 1;
+
+open(FH, "<", $ARGV[0]) || die "$!";
+
+while (<FH>) {
+	next unless /GEOMETRY/;
+	<FH>;
+
+	my @lines;
+
+	while (<FH>) {
+		last if /^$/;
+		push @lines, $_;
+	}
+
+	print scalar(@lines), "\n";
+	print "xyz", "\n";
+
+	foreach (@lines) {
+		$_ =~ s/A[0-9]+([a-zA-Z]+)[0-9]*/$1/;
+		print $_;
+	}
+}

.idea/modules.xml

@@ -227,11 +227,11 @@ Setting `enable_pbc` to `true` also sets `enable_cutoff` to `true`.
 
 ##### Periodic Box Size
 
-`periodic_box <x> <y> <z>`
+`periodic_box <x> <y> <z> <alpha> <beta> <gamma>`
 
-Default value: `30.0 30.0 30.0`
+Default value: `30.0 30.0 30.0 90.0 90.0 90.0`
 
-Unit: Angstrom
+Unit: Angstroms, degrees
 
 The smallest box dimension must be greater than `2 * swf_cutoff`.
 

.idea/vcs.xml

@@ -270,13 +270,16 @@ void print_matrix(size_t rows, size_t cols, const double *mat)
 	}
 }
 
-vec_t box_from_str(const char *str)
+six_t box_from_str(const char *str)
 {
-	vec_t box;
+	six_t box;
 
-	if (sscanf(str, "%lf %lf %lf", &box.x, &box.y, &box.z) < 3)
+	if (sscanf(str, "%lf %lf %lf %lf %lf %lf ", &box.x, &box.y, &box.z, &box.a, &box.b, &box.c) < 3)
 		error("incorrect box format");
 
-	vec_scale(&box, 1.0 / BOHR_RADIUS);
+	box.x *= 1.0 / BOHR_RADIUS;
+	box.y *= 1.0 / BOHR_RADIUS;
+	box.z *= 1.0 / BOHR_RADIUS;
+	// vec_scale(&box, 1.0 / BOHR_RADIUS);
 	return box;
 }

bin/cubegen.pl

@@ -117,7 +117,7 @@ void print_matrix(size_t, size_t, const double *);
 void check_fail(enum efp_result);
 void compute_energy(struct state *, bool);
 struct sys *parse_input(struct cfg *, const char *);
-vec_t box_from_str(const char *);
+six_t box_from_str(const char *);
 int efp_strcasecmp(const char *, const char *);
 int efp_strncasecmp(const char *, const char *, size_t);