#!/usr/bin/perl
#
# if your date command is different, change this too.
#
use lib '/home/ylin/Tools/perl/lib64/perl5/site_perl/5.8.5/x86_64-linux-thread-multi';
use Math::Cephes::Matrix qw(mat);
# 'mat' is a shortcut for Math::Cephes::Matrix->new
# perl script made by You Lin,
# Department of Physics,
# University of South Florida
# Dec, 2005
#Usage car2pwscf.pl file.car <option> <xmultu> <ymulti> <zmulti>
# <option>: : ATOMIC_POSITIONS in alat
# a: ATOMIC_POSITIONS in angstrom
# c: ATOMIC_POSITIONS in cell
# b: ATOMIC_POSITIONS in bohr
# d: Debug mode, shows inverse matrix and length of reci. vectors
# <xmultu> <ymulti> <zmulti>: Multiplications in x y z
$dateCmd = '/bin/date';
$PI = 3.14159265358979323846264;
# 4. * atan(1.);
$DTPI = $PI/180.;
$A0=0.5291772108; # Bohr radius in Angstrom
#-----------------------------------------------------------------------------
use constant pi => 3.14159265358979;
use constant EE => 2.718281828459045; # exp(1). Not 'e' because e is charge of electron
use constant Msun => 1.9891e+33;
use constant Lsun => 3.826e+33;
use constant Rsun => 6.9599e+10;
use constant c => 2.9979245620e+10; # speed of light
use constant G => 6.67259e-08; # gravitational constant
use constant e => 4.8032068e-10; # charge of electron
use constant h => 6.6260755e-27; # erg*s Planck's constant
use constant me => 9.1093897e-28; # mass of electron
use constant mp => 1.6726231e-24; # mass of proton
use constant alpha => 1/(h*c/(2*pi*e**2)); # fine structure constant
use constant sigmaT=> 6.6524616e-25; # Thomson cross section=8pi/3*re^2
use constant k => 1.380658e-16; # (erg/K) Boltzman constant
use constant NA => 6.0221367e+23; # mol^-1 Avogadro constant
use constant sigma => 5.67051e-5; # Stefan-Boltzmann constant
use constant keV => 1.602192e-9; use constant eV => 1.602192e-12; # erg
use constant pc => 3.085678e+18; use constant kpc => 1000*pc; use constant Mpc => 1000*kpc;
#-----------------------------------------------------------------------------
#$datafile = "@ARGV";
$datafile = $ARGV[0];
$opatm = $ARGV[1];
$XM = $ARGV[2];
$YM = $ARGV[3];
$ZM = $ARGV[4];
if ($ZM == 0) {
printf STDERR ("#WAR: No multiplication defined!\n");
printf STDERR ("#Default to 1 1 1\n");
$XM = 1 ; $YM = 1; $ZM = 1;
}
chop ($DateStamp = `$dateCmd +"%Y%m%d"`);
#$rec = join ( "\,", $DateStamp, $name );
open (DAT,"$datafile") || die "File $datafile Not found!\n";
$NUMATMS=0;
$NUMTYPES=0;
@ATMNAMES=("Empty");
while ($rec = <DAT>) {
@rec = split (" ",$rec);
if ($rec[0] =~ m/!|(PBC=)|(end)/i||$rec[0] eq "") {
} elsif ($rec[0] =~ m/PBC/i) {
@lattice=@rec;
} else {
$rec[0] =~ s/\d//g;
$NUMATMS ++;
$NUMTYPES=@ATMNAMES-1;
$FLAG=0;
for ($i=0;$i<$NUMTYPES;$i++) {
if ($rec[0] =~ $ATMNAMES[$i+1]) {
$FLAG=1;
break;
}
}
if ($FLAG == 0) {
push(@ATMNAMES,$rec[0]);
}
}
}
$NUMTYPES=@ATMNAMES-1;
if ($NUMTYPES == 0) {
printf STDERR ("ERR: NO ATOMS FOUND!\n");
exit;
}
$by=$lattice[2]*sin($DTPI*$lattice[6]);
$bx=$lattice[2]*cos($DTPI*$lattice[6]);
$cx=$lattice[3]*cos($DTPI*$lattice[5]);
$cy=($lattice[3]*$lattice[2]*cos($DTPI*$lattice[4])-$cx*$bx)/$by;
$cz=sqrt($lattice[3]*$lattice[3]-$cx*$cx-$cy*$cy);
my $LAT=$lattice[1];
my $M = mat([ [$LAT, $bx, $cx], [0, $by, $cy], [0, 0, $cz]]);
my $Mc = $M->coef;
my $I = $M->inv();
printf (" &system\n");
printf (" ibrav=0, celldm(1) =\$celldm, nat=%d, ntyp=$NUMTYPES,\n\n", $NUMATMS*$XM*$YM*$ZM);
if ($opatm eq "a") {
printf ("ATOMIC_POSITIONS {angstrom}\n");
} elsif ($opatm eq "c") {
printf ("ATOMIC_POSITIONS {crystal}\n");
} elsif ($opatm eq "b") {
printf ("ATOMIC_POSITIONS {bohr}\n");
} else {
printf ("ATOMIC_POSITIONS\n");
}
seek(DAT,0,0);
while ($rec = <DAT>) {
@rec = split (" ",$rec);
if ($rec[0] =~ m/!|(PBC=)|(end)/i||$rec[0] eq "") {
#printf ("%s\n",@rec);
} elsif ($rec[0] =~ m/PBC/i) {
#printf ("%s %20.9lf %20.9lf %20.9lf\n","Lattice",$rec[1],$rec[2],$rec[3]);
} else {
$rec[0] =~ s/\d//g;
for (my $i=0; $i<$XM; $i++) {
for (my $j=0; $j<$YM; $j++) {
for (my $k=0; $k<$ZM; $k++) {
my $T = [$rec[1],$rec[2],$rec[3]];
$R = [$T->[0]+$i*$Mc->[0][0]+$j*$Mc->[0][1]+$k*$Mc->[0][2],
$T->[1]+$i*$Mc->[1][0]+$j*$Mc->[1][1]+$k*$Mc->[1][2],
$T->[2]+$i*$Mc->[2][0]+$j*$Mc->[2][1]+$k*$Mc->[2][2]];
if ($opatm eq "a") {
} elsif ($opatm eq "c") {
$R = $I->mul($R);
$R = [$R->[0]/$XM, $R->[1]/$YM,$R->[2]/$ZM];
} elsif ($opatm eq "b") {
$R = [$R->[0]/$A0,$R->[1]/$A0,$R->[2]/$A0];
} else {
$R = [$R->[0]/$LAT,$R->[1]/$LAT,$R->[2]/$LAT];
}
printf ("%3s %20.9lf %20.9lf %20.9lf\n",$rec[0],$R->[0],$R->[1],$R->[2]);
}
}
}
}
}
close (DAT);
printf ("CELL_PARAMETERS\n");
for (my $i=0; $i<3; $i++) {
if ($i == 0) {$XYZM = $XM;}
elsif ($i == 1) {$XYZM = $YM;}
else {$XYZM = $ZM;}
printf (" %20.10lf %20.10lf %20.10lf\n",$Mc->[0][$i]/$LAT*$XYZM,
$Mc->[1][$i]/$LAT*$XYZM,$Mc->[2][$i]/$LAT*$XYZM);
}
printf ("\n Suggested \$celldm = %20.10lf Bohr.\n",$LAT/$A0);
if ($opatm eq "d") {
printf ("\n Inverse matrix:\n");
my $Ic = $I->coef;
for (my $i=0; $i<3; $i++) {
printf (" %20.10lf %20.10lf %20.10lf\n",$Ic->[0][$i]*$LAT/$XM,$Ic->[1][$i]*$LAT/$YM,$Ic->[2][$i]*$LAT/$ZM);
}
printf (" Reciprocal lattice length(1/Angstrom):\n");
for (my $i=0; $i<3; $i++) {
my $veclen=0.;
if ($i == 0) {$XYZM = $XM;}
elsif ($i == 1) {$XYZM = $YM;}
else {$XYZM = $ZM;}
for (my $j=0; $j<3; $j++) {
$veclen+=($Ic->[$i][$j]/$XYZM)*($Ic->[$i][$j]/$XYZM);
}
printf (" %20.15lf",sqrt($veclen));
}
printf ("\n");
}
|
#!/usr/bin/perl
#
# if your date command is different, change this too.
#
use lib '/home/ylin/Tools/perl/lib64/perl5/site_perl/5.8.5/x86_64-linux-thread-multi';
use Math::Cephes::Matrix qw(mat);
use Math::Trig;
# 'mat' is a shortcut for Math::Cephes::Matrix->new
# perl script made by You Lin,
# Department of Physics,
# University of South Florida
# Dec, 2005
#Usage pwscf2car.pl file.pwscf <latticeconstant(bohr)> <option>
# d: Debug mode
$dateCmd = '/bin/date';
$PI = 3.14159265358979323846264;
# 4. * atan(1.);
$DTPI = $PI/180.;
$A0=0.5291772108; # Bohr radius in Angstrom
#-----------------------------------------------------------------------------
use constant pi => 3.14159265358979;
use constant EE => 2.718281828459045; # exp(1). Not 'e' because e is charge of electron
use constant Msun => 1.9891e+33;
use constant Lsun => 3.826e+33;
use constant Rsun => 6.9599e+10;
use constant c => 2.9979245620e+10; # speed of light
use constant G => 6.67259e-08; # gravitational constant
use constant e => 4.8032068e-10; # charge of electron
use constant h => 6.6260755e-27; # erg*s Planck's constant
use constant me => 9.1093897e-28; # mass of electron
use constant mp => 1.6726231e-24; # mass of proton
use constant alpha => 1/(h*c/(2*pi*e**2)); # fine structure constant
use constant sigmaT=> 6.6524616e-25; # Thomson cross section=8pi/3*re^2
use constant k => 1.380658e-16; # (erg/K) Boltzman constant
use constant NA => 6.0221367e+23; # mol^-1 Avogadro constant
use constant sigma => 5.67051e-5; # Stefan-Boltzmann constant
use constant keV => 1.602192e-9; use constant eV => 1.602192e-12; # erg
use constant pc => 3.085678e+18; use constant kpc => 1000*pc; use constant Mpc => 1000*kpc;
#-----------------------------------------------------------------------------
#$datafile = "@ARGV";
$datafile = $ARGV[0];
$celldm = $ARGV[1]; # unit a.u.
if ($celldm == 0) {
printf STDERR ("ERR: No Lattice constant defined!\n");
exit;
}
$opatm="";
#chop ($DateStamp = `$dateCmd +"%Y%m%d"`);
chop ($DateStamp = `$dateCmd`);
#$rec = join ( "\,", $DateStamp, $name );
open (DAT,"$datafile") || die "File $datafile Not found!\n";
$NUMATMS=0;
$NUMTYPES=0;
@ATMNAMES=("Empty");
$FLAG=0;
$LINEINBK=0;
while ($rec = <DAT>) {
@rec = split (" ",$rec);
if ($rec[0] =~ m/&system/i ) {
$FLAG=1;
$LINEINBK=0;
} elsif ($rec[0] =~ m/ATOMIC_POSITIONS/i) {
$FLAG=2;
$opatm = $rec[1];
$LINEINBK=0;
} elsif ($rec[0] =~ m/CELL_PARAMETERS/i) {
$FLAG=3;
$LINEINBK=0;
} elsif ($rec[0] =~ m/K_POINTS/i||$rec[0] eq "") {
$FLAG=0;
$LINEINBK=0;
} else {
if ($FLAG == 3) {
if ($LINEINBK == 0) {
@latticex=@rec;
$LINEINBK ++ ;
} elsif ($LINEINBK == 1) {
@latticey=@rec;
$LINEINBK ++ ;
} elsif ($LINEINBK == 2) {
@latticez=@rec;
$LINEINBK ++ ;
}
} elsif ($FLAG == 2) {
$rec[0] =~ s/\d//g;
$NUMATMS ++;
$NUMTYPES=@ATMNAMES-1;
$FLAGIN=0;
for ($i=0;$i<$NUMTYPES;$i++) {
if ($rec[0] =~ $ATMNAMES[$i+1]) {
$FLAGIN=1;
break;
}
}
if ($FLAGIN == 0) {
push(@ATMNAMES,$rec[0]);
}
}
}
}
$NUMTYPES=@ATMNAMES-1;
if ($NUMTYPES == 0) {
printf STDERR ("ERR: NO ATOMS FOUND!\n");
exit;
}
my $M = mat([ [$latticex[0], $latticey[0], $latticez[0]],
[$latticex[1], $latticey[1], $latticez[1]],
[$latticex[2], $latticey[2], $latticez[2]]]);
my $Ind = mat([ [$celldm*$A0, 0, 0], [0, $celldm*$A0, 0], [0, 0, $celldm*$A0] ]);
$M = $M->mul($Ind);
my $Mc = $M->coef;
my $I = $M->inv();
my $Ic = $I->coef;
$Axx = $Ic -> [0][0]; $Axy = $Ic -> [0][1]; $Axz = $Ic -> [0][2];
$Ayx = $Ic -> [1][0]; $Ayy = $Ic -> [1][1]; $Ayz = $Ic -> [1][2];
$Azx = $Ic -> [2][0]; $Azy = $Ic -> [2][1]; $Azz = $Ic -> [2][2];
$a=sqrt($Azx*$Azx+$Azy*$Azy+$Azz*$Azz);
$Azx /= $a; $Azy /= $a; $Azz /= $a;
$a=sqrt($Ayx*$Ayx+$Ayy*$Ayy+$Ayz*$Ayz);
$Ayx /= $a; $Ayy /= $a; $Ayz /= $a;
$a=$Ayx*$Azx+$Ayy*$Azy+$Ayz*$Azz;
$Ayx=$Ayx-$a*$Azx; $Ayy=$Ayy-$a*$Azy; $Ayz=$Ayz-$a*$Azz;
$a=sqrt($Ayx*$Ayx+$Ayy*$Ayy+$Ayz*$Ayz);
$Ayx /= $a; $Ayy /= $a; $Ayz /= $a;
$a=$Axx*$Azx+$Axy*$Azy+$Axz*$Azz;
$b=$Axx*$Ayx+$Axy*$Ayy+$Axz*$Ayz;
$Axx=$Axx-$a*$Azx; $Axy=$Axy-$a*$Azy; $Axz=$Axz-$a*$Azz;
$Axx=$Axx-$b*$Ayx; $Axy=$Axy-$b*$Ayy; $Axz=$Axz-$b*$Ayz;
$a=sqrt($Axx*$Axx+$Axy*$Axy+$Axz*$Axz);
$Axx /= $a; $Axy /= $a; $Axz /= $a;
$A = mat([ [$Axx,$Axy,$Axz], [$Ayx,$Ayy,$Ayz],[$Azx,$Azy,$Azz]]);
my $NEWH = $A ->mul($M);
my $NEWHc = $NEWH -> coef;
$a=sqrt($NEWHc->[0][0]*$NEWHc->[0][0]);
$b=sqrt($NEWHc->[0][1]*$NEWHc->[0][1]+$NEWHc->[1][1]*$NEWHc->[1][1]);
$c=sqrt($NEWHc->[0][2]*$NEWHc->[0][2]+$NEWHc->[1][2]*$NEWHc->[1][2]+$NEWHc->[2][2]*$NEWHc->[2][2]);
$LAT=$a;
$gamma=$NEWHc->[0][1]/$b; $gamma=acos($gamma)/$DTPI;
$betta=$NEWHc->[0][2]/$c; $betta=acos($betta)/$DTPI;
$alpha=$NEWHc->[0][2]*$NEWHc->[0][1]+$NEWHc->[1][2]*$NEWHc->[1][1];
$alpha/=$b*$c; $alpha=acos($alpha)/$DTPI;
printf ("!BIOSYM archive 3\nPBC=ON\n!\n!DATE %s\n",$DateStamp);
printf ("PBC %9.4lf %9.4lf %9.4lf %9.4lf %9.4lf %9.4lf (P1)\n",
$a,$b,$c, $alpha, $betta, $gamma);
seek(DAT,0,0);
$FLAG=0;
$LINEINBK=0;
while ($rec = <DAT>) {
@rec = split (" ",$rec);
if ($rec[0] =~ m/&system/i ) {
$FLAG=1;
$LINEINBK=0;
} elsif ($rec[0] =~ m/ATOMIC_POSITIONS/i) {
$FLAG=2;
$LINEINBK=0;
} elsif ($rec[0] =~ m/CELL_PARAMETERS/i) {
$FLAG=3;
$LINEINBK=0;
} elsif ($rec[0] =~ m/K_POINTS/i||$rec[0] eq "") {
$FLAG=0;
$LINEINBK=0;
} else {
if ($FLAG == 3) {
if ($LINEINBK == 0) {
@latticex=@rec;
$LINEINBK ++ ;
} elsif ($LINEINBK == 1) {
@latticey=@rec;
$LINEINBK ++ ;
} elsif ($LINEINBK == 2) {
@latticez=@rec;
$LINEINBK ++ ;
}
} elsif ($FLAG == 2) {
$rec[0] =~ s/\d//g;
if ($opatm =~ m/crystal/) {
$S = [$rec[1],$rec[2],$rec[3]];
$R = $NEWH->mul($S);
} elsif ($opatm =~ m/angstrom/) {
$S = [$rec[1],$rec[2],$rec[3]];
$R = $A->mul($S);
} elsif ($opatm =~ m/bohr/) {
$S = [$rec[1]*$A0,$rec[2]*$A0,$rec[3]*$A0];
$R = $A->mul($S);
} elsif ($opatm =~ m/alat/) {
$S = [$rec[1]*$celldm*$A0,$rec[2]*$celldm*$A0,$rec[3]*$celldm*$A0];
$R = $A->mul($S);
} elsif ($opatm eq "") {
$S = [$rec[1]*$celldm*$A0,$rec[2]*$celldm*$A0,$rec[3]*$celldm*$A0];
$R = $A->mul($S);
}
$LINEINBK++;
$atmname = $rec[0]. "$LINEINBK";
printf ("%-5s %14.9lf %14.9lf %14.9lf XXXX 1 xx %-2s 0.000\n",
$atmname,$R->[0],$R->[1],$R->[2],$rec[0]);
}
}
}
printf ("end\nend\n");
close (DAT);
if ($ARGV[2] eq "d") {
printf ("\nCELL_PARAMETERS\n");
for (my $i=0; $i<3; $i++) {
printf (" %20.10lf %20.10lf %20.10lf\n",$NEWHc->[0][$i]/$LAT,$NEWHc->[1][$i]/$LAT,$NEWHc->[2][$i]/$LAT);
}
printf ("\n Lattice parameter = %20.10lf Ang.\n",$LAT);
}
exit;
|