Initial import of htcd system 1.0

Signed-off-by: Henrik Rydberg <rydberg@euromail.se>

1 files changed, 187 insertions, 0 deletions

diff --git a/src/labat/atomlab.cc b/src/labat/atomlab.cc
new file mode 100644
index 0000000..504a0d9
--- /dev/null
+++ b/src/labat/atomlab.cc
@@ -0,0 +1,187 @@
+/*************************************************************************
+ *
+ * HTCd - Copyright (C) 1998-2006 Henrik Rydberg
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include <mt/config.h>
+#include <fstream>
+#include <iostream>
+#include <unistd.h>
+#include <fcntl.h>
+#include <signal.h>
+#include <sys/stat.h>
+using namespace std;
+
+const mstring LABAT=htchome()+sref("/bin/labat");
+
+const int SLOTS=8;
+const int RUNLIMIT=600;
+static int SLOT;
+
+extern "C" {
+static void FreeSlot(int dummy)
+{
+  static int ok=0;
+  static char buf[256];
+  if(!ok) {
+    ok=1;
+    sprintf(buf,"rm -rf /tmp/labat%d/",SLOT);
+    system(buf);
+  }
+}
+}
+
+static int AllocSlot()
+{
+  struct stat FS;
+  char buf[256];
+  for(SLOT=0;SLOT<SLOTS;SLOT++) {
+    sprintf(buf,"/tmp/labat%d/",SLOT);
+    if(stat(buf,&FS)!=0) {
+      signal(SIGTERM,FreeSlot);
+      signal(SIGINT,FreeSlot);
+      signal(SIGPIPE,FreeSlot);
+      sprintf(buf,"mkdir /tmp/labat%d/",SLOT);
+      system(buf);
+      return 1;
+    }
+  }
+  return 0;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+struct CTLBasis {
+  int n,l,occ;
+  float e;
+  void SetEnergy(int Z) {
+    e=-0.5*float(Z*Z)/float(0.1+n*n*n+l*(l+1)*(l+2));
+  }
+};
+
+struct CTLData {
+  int atomz,its,basisn,exc;
+  float mix;
+  CTLBasis left[24],right[24];
+};
+
+static ostream& operator<<(ostream& out,const CTLData& ctl)
+{
+  out<<ctl.atomz<<" "<<0<<"\n";
+  out<<10<<" "<<40<<" "<<0.2<<"\n"; // Hermann-Skillman mesh
+  out<<ctl.mix<<" 0.00001 0.005\n"; // mix,error,?
+  out<<ctl.its<<" "<<0<<" "<<1<<" "<<0<<"\n"; // its ? step ?
+  out<<ctl.basisn<<" "<<2<<" "<<ctl.exc<<"\n"; // exc (C-A, B-H, GGA)
+  for(int i=0;i<ctl.basisn;i++) {
+    out<<ctl.left[i].n<<ctl.left[i].l<<0<<" "<<ctl.left[i].e
+       <<" "<<ctl.left[i].occ<<" "<<1<<" "<<0<<"\n"; 
+  }
+  for(int i=0;i<ctl.basisn;i++) {
+    out<<ctl.right[i].n<<ctl.right[i].l<<0<<" "<<ctl.right[i].e
+       <<" "<<ctl.right[i].occ<<" "<<1<<" "<<0<<"\n"; 
+  }
+  return out;
+}
+
+static void SetupLabat(mstring& resfile,int argc,char* argv[])
+{
+  CTLData ctl;
+  ctl.atomz=atoi(argv[1]);
+  ctl.its=atoi(argv[2]);
+  ctl.exc=atoi(argv[3]);
+  ctl.mix=atof(argv[4]);
+  resfile=argv[5];
+  ctl.basisn=0;
+  int nl,nr,elecs=0;
+  for(int i=0;i<6;i++) {
+    for(int k=0;k<3;k++) {
+      nl=int(argv[6+i][2*k]-'0');
+      nr=int(argv[6+i][2*k+1]-'0');
+      if(nl>0||nr>0) {
+        ctl.left[ctl.basisn].n=i+1;
+        ctl.left[ctl.basisn].l=k;
+        ctl.left[ctl.basisn].occ=nl;
+        ctl.left[ctl.basisn].SetEnergy(ctl.atomz);
+        ctl.right[ctl.basisn].n=i+1;
+        ctl.right[ctl.basisn].l=k;
+        ctl.right[ctl.basisn].occ=nr;
+        ctl.right[ctl.basisn].SetEnergy(ctl.atomz);
+        ctl.basisn++;
+      }
+      elecs+=nl+nr;
+    }
+  }
+  ctl.atomz+=elecs;
+  char buf[1024];
+  sprintf(buf,"/tmp/labat%d/atomctrl.dat",SLOT);
+  ofstream os(buf);
+  os<<ctl;
+  cout<<"Atomnumber: "<<ctl.atomz<<" Electrons: "<<elecs<<"<BR>\n";
+  cout<<"NOTE: if the program does >="<<ctl.its
+      <<" iterations, you must increase the number\n";
+  cout<<"of iterations and re-run, to assure convergence.<BR>\n";
+  cout<<"Also note that there is currently a maximum time limit of "<<
+    RUNLIMIT<<" seconds.<BR>\n";
+}
+
+static void RunLabat()
+{
+  char buf[1024];
+  sprintf(buf,"cd /tmp/labat%d/\nulimit -t %d\n%s",SLOT,RUNLIMIT,LABAT.c_str());
+  system(buf);
+}
+
+static void SaveLabat(const mstring& resfile)
+{
+  const char* docroot=getenv("DOCROOT");
+  if(docroot&&resfile.size()) {
+    char buf[1024];
+    sprintf(buf,"cp -rfp /tmp/labat%d/atomdens.dat %s/%s",
+            SLOT,docroot,resfile.c_str());
+    system(buf);
+    cout<<"<p>\n"
+        <<"Here you may download the density profile.\n"
+        <<"Example: Save the file as 'dens.dat', and load it from\n"
+        <<"matlab with 'load dens.dat'.\n"
+        <<"Then you can look at it with 'plot(dens(:,1),dens(:,2))'.\n"
+        <<"<p>\n"<<"<a href="<<resfile.c_str()<<">Get atomic density</a>\n";
+  }
+}
+
+main(int argc,char* argv[])
+{
+  if(argc<12) {
+    cerr<<"Usage: "<<argv[0]<<" <ion> <its> <exc> <mix> <resfile> <ssppdd>x6\n";
+    exit(-1);
+  }
+  mstring resfile;
+  char buf[1024];
+  if(AllocSlot()) {
+    cout<<"<h4>Running</h4>\n"
+        <<"<PRE>\n";
+    cout.flush();
+    SetupLabat(resfile,argc,argv);
+    RunLabat();
+    cout<<"</PRE>\n"
+        <<"<h4>Done</h4>\n";
+    SaveLabat(resfile);
+    FreeSlot(0);
+  }
+  else {
+    cout<<"No vacant slots at the moment. Try again later.\n";  
+  }
+}