Commit 11c86f2f authored by Alberto Nocera's avatar Alberto Nocera

added script for Cheby, only chains

parent c08baba7
#!/usr/bin/perl
use strict;
use warnings;
use Math::Trig;
use lib ".";
use OmegaUtils;
use Getopt::Long qw(:config no_ignore_case);
my $usage = "-f Input [-M mForQ] [-S site] [-p] [-r] [-z]\n";
my ($Input,$site,$m,$GlobalNumberOfSites);
my ($siteForSpectrum,$mForQ,$isPeriodic,$mMax,$wantsRealPart);
my $zeroAtCenter = 0;
GetOptions('f=s' => \$Input,
'S:i' => \$siteForSpectrum,
'm:i' => \$mForQ,
'p' => \$isPeriodic,
'M:i' => \$mMax,
'r' => \$wantsRealPart,
'z' => \$zeroAtCenter) or die "$usage\n";
(defined($Input) and defined($isPeriodic)) or die "$0: USAGE: $usage\n";
my ($omega0,$total,$omegaStep,$centralSite,$Wstar,$epsilont,$JacksOrLorentz);
my $geometryName;
my $geometryLeg = 1;
my $orbitals = 1;
my $hptr = {"#OmegaBegin" => \$omega0,
"#OmegaTotal" => \$total,
"#OmegaStep" => \$omegaStep,
"GeometryKind" => \$geometryName,
"LadderLeg" => \$geometryLeg,
"Orbitals" => \$orbitals,
"ChebyshevWstar" => \$Wstar,
"ChebyshevEpsilon" => \$epsilont,
"TSPSites 1" => \$centralSite,
"#JacksonOrLorentz" => \$JacksOrLorentz,
"TotalNumberOfSites" => \$GlobalNumberOfSites};
OmegaUtils::getLabels($hptr,$Input);
$hptr->{"isPeriodic"} = $isPeriodic;
$hptr->{"mMax"} = $mMax;
$hptr->{"centralSite"} = $centralSite;
my $L = $GlobalNumberOfSites;
my $time = 0.0;
my @vM;
my $siteC = $centralSite;
$Input =~s/.inp//;
open(FIN,"runFor$Input.cout") or die "$0: Cannot open runForinput3511.cout : $!\n";
print STDERR "runFor$Input.cout opened\n";
my $status;
while (<FIN>) {
#next unless (/ALL OPERATORS HAVE BEEN APPLIED/);
if (/P1/ and /gs/) {
$status="p1";
} elsif (/P2/ and /gs/) {
$status="p2";
} elsif (/P2/ and /gs/ and /PsiApp:/) {
next;
} else {
next;
}
chomp;
my @temp = split;
if ($temp[0] eq 'PsiApp:') {next;}
my $p = $temp[0];
$time = $temp[2];
#print STDERR "$temp[0] $temp[1] $temp[2] $temp[3] $temp[4]\n";
$vM[$p+$L*$siteC+$L*$L*$time] = $temp[1] if ($status eq "p1");
}
close(FIN);
print STDERR "MaxTime=$time\n";
for (my $p = 1; $p<$L-1; $p++) {
for (my $i =0; $i<$time; $i++) {
#print STDERR " $vM[$p+$L*$siteC+$L*$L*$i]\n";
}
}
my @dampG;
$dampG[0] = 1.0;
if ($JacksOrLorentz==1) {
for (my $n = 1; $n<$time; $n++) {
my $num = ($time - $n + 1.0) * cos(pi*$n/($time+1.0)) + sin(pi*$n/($time+1.0)) * cot(pi/($time+1.0));
my $den = ($time + 1.0);
$dampG[$n] = $num/$den;
}
} else {
for (my $n = 1; $n<$time; $n++) {
my $num = sinh(4.0*(1.0-$n/$time));
my $den = sinh(4.0);
$dampG[$n] = $num/$den;
}
}
#print STDERR "Damping Factor\n";
#for (my $i =0; $i<$time; $i++) {
#print STDERR "$i\t$dampG[$i]\n";
#}
# Preparing the Chebyshev Polynomials
my $Wprime = 1.0-0.5*$epsilont;
my $scalea = $Wstar/(2.0*$Wprime);
my @factor;
my @ChebyPolyXfactor;
for (my $om =0; $om<$total; $om++) {
for (my $n = 1; $n<$time; $n++) {
$ChebyPolyXfactor[$om+$total*$n] = 0.0;
}
}
for (my $om =0; $om<$total; $om++) {
my $omega = $omega0+$om*$omegaStep;
# Scaling
my $omegaPrime = ($omega/$scalea)-$Wprime;
$factor[$om] = 1.0/(pi*sqrt(1.0-$omegaPrime*$omegaPrime));
for (my $n = 1; $n<$time; $n++) {
$ChebyPolyXfactor[$om+$total*$n] = $factor[$om]*cos($n*acos($omegaPrime));
}
}
# Constructing the Real Space Structure Factor
my @vMdamped;
for (my $om =0; $om<$total; $om++) {
for (my $p = 0; $p<$L; $p++) {
$vMdamped[$p+$L*$siteC+$L*$L*$om] = 0.0;
}
}
for (my $om =0; $om<$total; $om++) {
for (my $p = 1; $p<$L-1; $p++) {
# adding the zeroth Cheby momentum $n==0
$vMdamped[$p+$L*$siteC+$L*$L*$om] += $factor[$om]*$vM[$p+$L*$siteC];
# adding all the others
for (my $n = 1; $n<$time; $n++) {
$vMdamped[$p+$L*$siteC+$L*$L*$om] += 2.0*$dampG[$n]*$ChebyPolyXfactor[$om+$total*$n]*$vM[$p+$L*$siteC+$L*$L*$n];
}
}
}
if (!defined($mForQ)) {$mForQ = $L;}
my $outSpectrum = "outSpectrum.L=$L.gnuplot";
open(FOUTSPECTRUM,"> $outSpectrum") or die "$0: Cannot write to $outSpectrum : $!\n";
for (my $om =0; $om<$total; $om++) {
my $omega = $omega0+$om*$omegaStep;
for (my $nk =0; $nk<$mForQ; $nk++) {
my $sum1 = 0;
my $k = 2.0*($nk)*pi/($mForQ);
for (my $p = 0; $p<$L; $p++) {
$sum1 += (1/$L) * cos($k*($p-$siteC)) * $vMdamped[$p+$L*$siteC+$L*$L*$om];
}
print FOUTSPECTRUM "$k $omega $sum1\n";
}
#print FOUTSPECTRUM "\n";
}
close(FOUTSPECTRUM);
print STDERR "$0: Spectrum written to $outSpectrum\n";
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment