procOmegas: Grand Canonical Spectra for physical chain

scripts/OmegaUtils.pm

@@ -206,6 +206,9 @@ sub getGeometryDetails
 	if ($name eq "chain") {
 		$factor = 0.5;
 		die "$0: Chain does not have ky != 0\n" if (defined($my) and $my != 0)
+	} elsif ($subname eq "GrandCanonical" and $name eq "ladder") {
+		$factor = 0.5;
+		die "$0: Chain does not have ky != 0\n" if (defined($my) and $my != 0);
 	} elsif ($name eq "ladder" || $subname eq "average") {
 		$leg = $geometry->{"leg"};
 		$factor = 0.25;
@@ -254,6 +257,10 @@ sub fourier
 		return fourierChain($f, $v, $hptr);
 	}
 
+	if ($name eq "ladder" and $subname eq "GrandCanonical") {
+		return fourierChainGC($f, $v, $hptr);
+	}
+
 	if ($name eq "ladder") {
 		return fourierLadder($f, $v, $geometry->{"leg"}, $hptr);
 	}
@@ -575,6 +582,45 @@ sub honeyFourierFactor
 	return (cos($arg), sin($arg));
 }
 
+sub fourierChainGC
+{
+	my ($f, $v, $hptr) = @_;
+	my $n = scalar(@$v);
+
+	my $mMax = $hptr->{"mMax"};
+	my $isPeriodic = $hptr->{"isPeriodic"};
+	my $centralSite = $hptr->{"centralSite"};
+	my $nOver2 = int($n/2);
+
+	die "$0: FATAL: ChainGC central site is odd\n" if ($centralSite & 1);
+
+	if (!$isPeriodic) {
+		my $b = ($centralSite != $nOver2);
+		if ($b && $centralSite != $nOver2 - 2) {
+			die "$0: FATAL ChainGC: wrong central site $centralSite\n";
+		}
+	}
+	
+	my $cSite = int($centralSite/2);
+	my $numberOfQs = (defined($mMax)) ? $mMax : $nOver2;
+	for (my $m = 0; $m < $numberOfQs; ++$m) {
+		my @sum = (0,0);
+		my $q = getQ($m, $numberOfQs, $isPeriodic);
+		for (my $ii = 0; $ii < $n; $ii += 2) {
+			my $i = int($ii/2);
+			my $ptr = $v->[$ii];
+			my @temp = @$ptr;
+			my $arg = $q*($i - $cSite);
+			my $carg = cos($arg);
+			my $sarg = sin($q*($i + 1))*sin($q*($cSite + 1));
+			my $cOrSarg = ($isPeriodic) ? $carg : $sarg;
+			$sum[0] += $temp[0]*$cOrSarg;
+			$sum[1] += $temp[1]*$cOrSarg;
+		}
+
+		$f->[$m] = \@sum;
+	}
+}
 
 sub writeFourier
 {
@@ -582,7 +628,7 @@ sub writeFourier
 	my $subname = $geometry->{"subname"};
 	my $isPeriodic = $geometry->{"isPeriodic"};
 
-	if ($geometry->{"name"} eq "chain") {
+	if ($geometry->{"name"} eq "chain" || $subname eq "GrandCanonical") {
 		return writeFourierChain($array,$f, $isPeriodic);
 	}
 

scripts/procOmegas.pl

@@ -49,7 +49,7 @@ my $hptr = {"#OmegaBegin" => \$omega0,
             "#OmegaOffset" => \$omegaOffset,
             "#Lx" => \$lx,
             "#Ly" => \$ly,
-	    "#options" => \$options,
+            "#options" => \$options,
             "GeometryKind" => \$geometryName,
             "GeometrySubKind" => \$geometrySubName,
             "LadderLeg" => \$geometryLeg,
@@ -67,6 +67,7 @@ $centralSite = getCentralSite($templateInput, $isAinur);
 if ($isAinur) {
 	$geometryName =~ s/[\";]//g;
 	$GlobalNumberOfSites =~ s/;//g;
+	$geometrySubName =~ s/[\";]//g;
 }
 
 $hptr->{"isPeriodic"} = $isPeriodic;
@@ -197,6 +198,7 @@ sub procCommon
 
 	my @qValues;
 	OmegaUtils::fourier(\@qValues,\@spaceValues,$geometry,$hptr);
+	print LOGFILEOUT "$0: Number of k values ".scalar(@qValues)."\n";
 	OmegaUtils::writeFourier($array,\@qValues,$geometry);
 }
 
@@ -579,7 +581,7 @@ sub getRealOrImagData
 	my @temp;
 	my $n = scalar(@$d);
 	my $start = 1;
-	if ($geometry->{"name"} eq "ladder") {
+	if ($geometry->{"name"} eq "ladder" and $geometry->{"subname"} ne "GrandCanonical") {
 		my $leg = $geometry->{"leg"};
 		$n = int($n/$leg);
 		$start += $qyIndex*$n;