535. Briefly: ACS journal latex template -- achemso

There's no outright template for the preparation of tex files for ACS journals. However, there are two files provided: a bibtex styles file, achemso.bst, and a macro file, achemso.cls. Both are available in the texlive-latex-extra package on debian.

Luckily there's a lot of information online, including a demo file: ftp://ftp.dante.de/tex-archive/biblio/bibtex/contrib/achemso/

Anyway, here's a simple file that can act as a template. Have a look at the demo file at ftp.dante.de for a much more exhaustive example, including how to use schemes and insert references.

You can figure out the journal abbreviations from their URLs. Otherwise, page 4 in this pdf has a list: http://www.tug.org/texlive/Contents/live/texmf-dist/doc/latex/achemso/achemso.pdf

Anyway:

\documentclass[journal=inoraj, layout=twocolumn]{achemso}
\usepackage[latin1]{inputenc}
\usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\author{Vera Hill}
\affiliation{Department of Chemistry, This University, This Country}
\email{my@email}
\author{I. Lindqvist}
\affiliation{Department of Chemistry, Another University, That Country}
\title{A quick, non-exhaustive tex template}

\begin{document}
\maketitle

%\begin{tocentry} %graphical TOC
%\includegraphics{example.eps}
%TOC text goes here
%\end{tocentry}

\begin{abstract}
The ACS should provide a simple template. They don't, so I do.
\end{abstract}

\section{Introduction}
\section{Results and Discussion}

%\begin{figure}
% \includegraphics{graphic}
% \caption{A figure}
% \label{fig:example}
%\end{figure}


\subsection{References}
\section{Experimental}
\begin{acknowledgement}
VH thanks the internet. IL thanks the electron.
\end{acknowledgement}

\begin{suppinfo}
See supporting information for additional experimental details.
\end{suppinfo}

\end{document}

---

Example:

534. Adding new options to ECCE -- adding solvation models for G03 (G09)

To minimize the amount of manual editing of my input files I've started to modify the menus in ECCE. One thing that was missing before was the option of selecting a solvation model for Gaussian jobs. I've added that now.

I will eventually submit my changes upstreams, but for now I'll just post a clunky list of changes that I've made:
There are two files to edit:

apps/scripts/codereg/ged03theory.py.

Add the following right after the DFT options section and before the MP options section (ca line number 273 in my case)

#Theory options solvation -CAO

            sovSizer = EcceBoxSizer(self,
                                    label = "Solvation",
                                    cols = 2)
            sovLeftSizer = EcceVBoxSizer()
            sovRightSizer = EcceVBoxSizer()


           #Use solvation
            self.useSCRF = EcceCheckBox(self, #useCosmo
                                         label = " Use SCRF"
                                         name = "ES.Theory.SCF.UseSCRF",
                                         default = False)
            sovLeftSizer.AddWidget(self.useSCRF,
                                   border = EcceGlobals.BorderDefault)
           #SCRF type
            scrfChoice = ["PCM",
                             "CPCM",
                             "IPCM",
                             "SCIPCM",
                             "SMD",
                             "Dipole"]
            self.scrf = EcceComboBox(self,
                                        choices = scrfChoice,
                                        name = "ES.Theory.SCF.SCRF",
                                        label = "SCRF type:",
                                        default = 0)
            sovRightSizer.AddWidget(self.scrf,
                                    border = EcceGlobals.BorderDefault)
           #Solvent type

            solventChoice = ["Water",
                             "Acetonitrile",
                             "Methanol",
                             "Ethanol",
                             "Manual",
                             "Benzene",
                             "Chloroform",
                             "Diethylether",
                             "Dichloromethane",
                             "Dichloroethane",
                             "Carbontetrachloride",
                             "Toluene",
                             "Chlorobenzene",
                             "Nitromethane",
                             "Heptane",
                             "Aniline",
                             "Acetone",
                             "Tetrahydrofuran",
                             "Dimethylsulfoxide",

                             "Argon",
                             "Krypton",
                             "Xenon",
                             "n-Octanol",
                             "1-Butanol"
                             "Cyclohexane",
                             "Isoquinoline",
                             "Quinoline",
                             ]   
            self.solvent = EcceComboBox(self,
                                        choices = solventChoice,
                                        name = "ES.Theory.SCF.Solvent",
                                        label = "Solvent:",
                                        default = 0)
            sovRightSizer.AddWidget(self.solvent,
                                    border = EcceGlobals.BorderDefault)

            self.scrfDielec = EcceFloatInput(self,
                                              default = 78.4,
                                              name = "ES.Theory.SCF.Dielectric",
                                              label = "Dielectric Constant:",
                                              hardRange = "(0..)",
                                              unit = "Debye")
            sovRightSizer.AddWidget(self.scrfDielec,
                                    border = EcceGlobals.BorderDefault)

            sovSizer.AddWidget(sovLeftSizer,
                               flag = wx.ALL)
            sovSizer.AddWidget(sovRightSizer,
                               flag = wx.ALL)
            self.panelSizer.Add(sovSizer)
#End theory options solvation -CAO

Also add the following to the def CheckDependency(self) block, right before the "if (EcceGlobals.Category == "MP" or" line:

#       SCRF solvation -- CAO 
        if EcceGlobals.Category == "DFT" or EcceGlobals.Category == "SCF":   
            self.solvent.Enable(self.useSCRF.GetValue())
            self.scrf.Enable(self.useSCRF.GetValue())
            self.scrfDielec.Enable(self.useSCRF.GetValue() and 
                                    self.solvent.GetValue() == "Manual")
#       end SCRF

Edit apps/scripts/parsers/ai.gauss03 and put the following somewhere in the file.

##############################################################################
#
#  Description:
#      SCRF options field -CAO
#       
############################################################################## 
sub SCRFOptions {
  my($options,$result);

  $result = "";
  $options = "";

#scrf type
if ($AbiDict{"ES.Theory.SCF.UseSCRF"}) {
 if ($AbiDict{"ES.Theory.SCF.SCRF"} eq "" ) {
  $options .= "PCM,";
  }
 else {$options .= ($AbiDict{"ES.Theory.SCF.SCRF"}).","};

  
# solvent/dielectric 
 if ($AbiDict{"ES.Theory.SCF.Solvent"} eq "Manual") {
     if (defined($AbiDict{"ES.Theory.SCF.Dielectric"})) {
       $options .= "Dielectric=".$AbiDict{"ES.Theory.SCF.Dielectric"}."";
      } }
      elsif ($AbiDict{"ES.Theory.SCF.Solvent"} eq "") {
    $options .= "Solvent=water"; }
      else { $options .= "Solvent=".$AbiDict{"ES.Theory.SCF.Solvent"};
   }

  if ($options ne "") {
    $result =  "SCRF=(";
    $result .= $options;
    $result .= ") ";
  }

  return  $result;
}
}

Also add the following right below the "$route .= &SCFOptions;" line:

$route .= &SCRFOptions;

And you're done:

533. Adding a new Exchange/Correlation functional to ECCE. M06 for G03 (G09)

Nothing fancy here. I'd just like to be able to select the M06 functional, which is available in G09, using the ECCE interface.

At some point in the future I might set up separate files for G09, but for now I run G09 using the G03 files in ECCE.

See here and here for how to add basis sets to ECCE.

You'll need to edit 2-3 files. I did this in ECCE v7.0.

Edit apps/scripts/codereg/ged03theory.py

185             xcFuncChoice = ["None",
186                             "M06 (hybrid)",
187                             "SVWN 5 (local)",
188                             "SVWN 1/RPA (local)",
189                             "BLYP (nonlocal)",

Edit apps/scripts/parsers/ai.gauss03

635   if ($xcFun eq "SVWN 1/RPA (local)") {# convert our terminology to Gaussian's
 636     $result = "SVWN";
 637   } elsif ($xcFun eq "SVWN 5 (local)") {
 638     $result = "SVWN5";
 639   } elsif ($xcFun eq "M06 (hybrid)") {
 640     $result = "M06";
 641   } elsif ($xcFun eq "BLYP (nonlocal)") {

Editing apps/scripts/parsers/Gaussian-03.expt isn't necessary to get this to work, but I did it anyway. Not sure if it 'does' anything:

557       if ($method =~ /b.*3lyp/  ||
 558           $method =~ /m06/      ||
 559           $method =~ /bhandh/   ||
 560           $method =~ /vscx/     ||
 561           $method =~ /hcth407/  ||
 562           $method =~ /hcth147/  ||

You should now be able to select M06 in the Theory Details dialogue:

532. TEMPer temperature monitoring USB stick on Debian Wheezy

Because the air conditioning in my office has a habit of turning itself off, and since I'm running my beowulf cluster in there, and since it's Australia, I've become interested in monitoring the temperature in my office.

The USB stick itself looks nothing special, so here's the card from the box it came in.

A colleague of mine got a TEMPer thermometer USB (0c45:7401 Microdia) back when he didn't have any air conditioning at all in his office and wanted to prove to the university that the temperature got so high that it was impossible for him to do any work on some days. He's now got air conditioning.

Anyway, plugging in the USB stick got me the following:
* /dev/hidraw5 and /dev/hidraw6 get created

* DMESG shows

[441126.932728] usb 2-4.2: new low-speed USB device number 11 using ehci-pci
[441127.025790] usb 2-4.2: New USB device found, idVendor=0c45, idProduct=7401
[441127.025803] usb 2-4.2: New USB device strings: Mfr=1, Product=2, SerialNumber=0
[441127.025811] usb 2-4.2: Product: TEMPerV1.2
[441127.025818] usb 2-4.2: Manufacturer: RDing
[441127.030229] input: RDing TEMPerV1.2 as /devices/pci0000:00/0000:00:02.1/usb2/2-4/2-4.2/2-4.2:1.0/input/input24
[441127.030516] hid-generic 0003:0C45:7401.000F: input,hidraw5: USB HID v1.10 Keyboard [RDing TEMPerV1.2] on usb-0000:00:02.1-4.2/input0
[441127.033234] hid-generic 0003:0C45:7401.0010: hiddev0,hidraw6: USB HID v1.10 Device [RDing TEMPerV1.2] on usb-0000:00:02.1-4.2/input1

* lsusb shows

Bus 002 Device 011: ID 0c45:7401 Microdia 

Searching online for 0c45:7401 brought up this cheesily title post: http://www.linuxjournal.com/content/temper-pi

From that post:

 If instead dmesg says this:
[snip]
and lsusb says:
$ lsusb
Bus 001 Device 005: ID 0c45:7401 Microdia
then congratulations, you have the new TEMPer probe and will have to use completely different software. 

While that sounds as if you'll have continue searching for a new how-to, in fact the entire post is about that particular version. So, I followed the instructions at Linux Journal -- I'll just offer my step by step version of it here with some added detail:

sudo apt-get install python-usb python-setuptools snmpd git
sudo easy_install snmp-passpersist
mkdir ~/tmp
cd ~/tmp
git clone git://github.com/padelt/temper-python.git
cd temper-python/
sudo python setup.py install

At this point I could get a temperature reading by doing:

$ sudo temper-poll 

Found 1 devices
Device #0: 24.4°C 75.9°F

But running stuff as root is unsatisfying, so I created a UDEV rule:

$ sudo vim /etc/udev/rules.d/80-temper.rules

SUBSYSTEMS=="usb", ATTRS{idVendor}=="0c45", ATTRS{idProduct}=="7401", GROUP="users", MODE="0666"

I then unplugged the USB stick, did

sudo service udev restart

and plugged it back in.

$ temper-poll 
Found 1 devices
Device #0: 25.8°C 78.3°F

Sweet.
Finally, I set up a cronjob that would check the temperature, update a plot and put it in my Dropbox:

$ crontab -e

*/2 *   * * *   sh /home/me/temper.sh

where temper.sh looks like this:

temp=`/usr/local/bin/temper-poll |grep Device|gawk '{print $3}'|sed 's/°C//'`
when=`date +%s`
thetime=`date +%D' '%T`
if [ -n "$temp" ]; then
        echo $when $temp $thetime>> /home/me/temper.dat
fi

gnuplot /home/andy/temper.gplt
cp /home/me/temper.eps /home/me/Dropbox

The temper.gplt script looks like this:

set term postscript eps enhanced colour
set output 'temper.eps'

unset key
set ylabel 'Temperature (Celsius)'
set border 3
set xtics nomirror
set ytics nomirror
unset xlabel
set xdata time

set multiplot
        set size 0.5,0.45

        set origin 0,0.05
                set timefmt "%H:%M:%S"
                set title 'Daily'
                set xtics 30000
                plot 'temper.dat' u 4:2 w points pt 1 ps 0.15

        set origin 0.5,0.05
                set title 'By Day'
                set timefmt "%m/%d/%y"
                set xtics 100000#0
                plot 'temper.dat' u 3:2 w points pt 2 ps 0.5

        set size 1.0,0.5
        set origin 0.0,0.5
                set timefmt "%m/%d/%y %H:%M:%S"
                set xtics 30000
                set title 'Log'
                plot 'temper.dat' u 3:2 w lines

and the plot looks like this:

Temperature in a lab at a leading Australian research institute. In five years they have not been able to fix the air conditioning.  On 01/04 someone pushed a cardboard box against the sensor which lead to a slower change in temperature.

531. Briefly: NWChem 6.3 -- issues with planewave (PSPW) module and AMD FX8150 and 8350 CPUs

This is more of an announcement or warning than a proper blog post:

Both FX8350 and FX8150 have trouble running the pspw module causing the calculation to lead to exploding structures:
http://www.nwchem-sw.org/index.php/Special:AWCforum/st/id1059/Issue_with_pspw_using_nwchem_6.3....html

My other nodes have no trouble running the job in question. Also, the issue was only found in nwchem 6.3 -- nwchem 6.1.1 worked fine. So it's not an FX83x50 related fault per se.

Again, see the post at the nwchem-sw.org site for more information.

530. Briefly: Adding a new entry to Default applications in Gnome 3. Example using Firefox

I'm tired of google chrome/chromium -- for some reason more and more websites are rendering incorrectly in it. Part of the reason is because I refuse to allow just any website to set cookies, but that can't explain all instances (e.g. I'm having major issues with any elsevier journals). I'm also tired of google in general, especially after having spent some time with google plus.
Anyway, I recently showed how to install firefox: http://verahill.blogspot.com.au/2013/11/528-briefly-setting-up-64-bit-pre-built.html

I'll show here to set it as a selectable application in the System Settings/Details

At the beginning the following selections are available:

The key to adding a new 'Default Application' is simply making sure that it appears in the MIME file type associations. And one way to do that is to create a .desktop file and use update-desktop-database to read it:

sudo cp /usr/share/applications/iceweasel.desktop /usr/share/applications/firefox.desktop
sudo sed -i 's/Iceweasel/Firefox/g' /usr/share/applications/firefox.desktop
sudo sed -i 's/Exec=iceweasel/Exec=firefox25/g' /usr/share/applications/firefox.desktop
sudo sed -i 's/Icon=iceweasel/Icon=firefox/g' /usr/share/applications/firefox.desktop
sudo update-desktop-database

Once that's done we get the following:

 You can then set up a shortcut launcher, e.g.

(you could of course just have it execute the command directly, but what's the challenge in that?)

529. Briefly: Error Writing spool: NT_STATUS_DISK_FULL

I recently had trouble printing on a networked printer at work, where we use 'Papercut' to share printers -- basically you submit your job, give your credentials, then run over to a printer and release the job.

Anyway, I suddenly had issues printing:

The solution:
I'm not entirely sure what fixed it, but here's what I did

lpstat showed a number of jobs that had been submitted to the printer, but couldn't be released:

me@beryllium:~/Downloads$ lpstat
global-mfp-1166         me              591872   Mon 11 Nov 2013 12:36:29 EST
global-mfp-1167         me              993280   Mon 11 Nov 2013 12:36:44 EST
global-mfp-1168         me             2014208   Mon 11 Nov 2013 12:36:59 EST
global-mfp-1169         me              871424   Mon 11 Nov 2013 12:37:17 EST
global-mfp-1170         me              573440   Mon 11 Nov 2013 12:37:31 EST
global-mfp-1171         me             1199104   Mon 11 Nov 2013 12:37:51 EST
global-mfp-1172         me              183296   Mon 11 Nov 2013 12:38:02 EST
global-mfp-1173         me              491520   Mon 11 Nov 2013 12:38:19 EST
global-mfp-1174         me             2035712   Mon 11 Nov 2013 12:38:38 EST
global-mfp-1175         me             2035712   Mon 11 Nov 2013 12:39:54 EST
global-mfp-1176         me              635904   Mon 11 Nov 2013 12:41:54 EST
global-mfp-1177         me              148480   Mon 11 Nov 2013 16:29:58 EST

I preceded to cancel all the jobs:

me@beryllium:~/Downloads$ cancel global-mfp-1166
me@beryllium:~/Downloads$ lpstat
global-mfp-1167         me              993280   Mon 11 Nov 2013 12:36:44 EST
global-mfp-1168         me             2014208   Mon 11 Nov 2013 12:36:59 EST
global-mfp-1169         me              871424   Mon 11 Nov 2013 12:37:17 EST
global-mfp-1170         me              573440   Mon 11 Nov 2013 12:37:31 EST
global-mfp-1171         me             1199104   Mon 11 Nov 2013 12:37:51 EST
global-mfp-1172         me              183296   Mon 11 Nov 2013 12:38:02 EST
global-mfp-1173         me              491520   Mon 11 Nov 2013 12:38:19 EST
global-mfp-1174         me             2035712   Mon 11 Nov 2013 12:38:38 EST
global-mfp-1175         me             2035712   Mon 11 Nov 2013 12:39:54 EST
global-mfp-1176         me              635904   Mon 11 Nov 2013 12:41:54 EST
global-mfp-1177         me              148480   Mon 11 Nov 2013 16:29:58 EST
me@beryllium:~/Downloads$ cancel global-mfp-1167
me@beryllium:~/Downloads$ cancel global-mfp-1168
me@beryllium:~/Downloads$ cancel global-mfp-1169
me@beryllium:~/Downloads$ cancel global-mfp-1170
me@beryllium:~/Downloads$ cancel global-mfp-1171
me@beryllium:~/Downloads$ cancel global-mfp-1172
me@beryllium:~/Downloads$ cancel global-mfp-1173
me@beryllium:~/Downloads$ cancel global-mfp-1174
me@beryllium:~/Downloads$ cancel global-mfp-1175
me@beryllium:~/Downloads$ cancel global-mfp-1176
me@beryllium:~/Downloads$ cancel global-mfp-1177
me@beryllium:~/Downloads$ lpstat
me@beryllium:~/Downloads$ 

That didn't remove the error message, however. Opening Printers in Gnome(3) showed that the printer with an issues was set to 'off'. Unlocking and changing it to 'on' resolved the issue.

528. Briefly: Setting up 64 bit pre-built firefox (25) binaries on debian

/usr/local/lib is just a suggestion (and probably not a very good one).

Anyway, here's the quick and easy way to get set up with 64 bit firefox. Start it by calling firefox25

mkdir ~/tmp
cd ~/tmp
wget ftp://ftp.mozilla.org/pub/firefox/releases/25.0/linux-x86_64/en-GB/firefox-25.0.tar.bz2
tar xvf firefox-25.0.tar.bz2 
cd firefox/
sudo mkdir /usr/local/lib/firefox-25.0
sudo cp * -R /usr/local/lib/firefox-25.0
sudo ln -s /usr/local/lib/firefox-25.0/firefox /usr/bin/firefox25

527. Briefly: setting up thunderbird 24.1.0 on Debian (binaries)

Not much to say, other than that building thunderbird is a bit more complex these days than simply doing a configure/make/make install pass. For once I decided that rolling my own version wasn't worth it, and grabbed the pre-built binaries instead.

Here's a very brief description of how to get them set up:

mkdir ~/tmp
cd ~/tmp
wget ftp://ftp.mozilla.org/pub/mozilla.org/thunderbird/releases/24.1.0/linux-x86_64/en-GB/thunderbird-24.1.0.tar.bz2
tar xvf thunderbird-24.1.0.tar.bz2 
cd thunderbird/
sudo mkdir /usr/local/lib/thunderbird-24.1.0
sudo cp * -R /usr/local/lib/thunderbird-24.1.0
sudo rm /usr/local/bin/thunderbird
sudo ln -s /usr/local/lib/thunderbird-24.1.0/thunderbird /usr/local/bin/thunderbird

526. New release of NWChem 6.3 out (17th of October 2013)

There's recently a new release of nwchem 6.3 (release 2). As usual there's no public message, no release notes or anything that actually informs you as to whether there are critical bug fixes, new functionality or anything else.

The new version can be found here: http://www.nwchem-sw.org/download.php?f=Nwchem-6.3.revision2-src.2013-10-17.tar.gz

I'm not competent in telling you whether you should upgrade or not, but here's a list over the changed files:

nwchem-6.3.revision2-src.2013-10-17/INSTALL
nwchem-6.3.revision2-src.2013-10-17/src/config/makefile.h
nwchem-6.3.revision2-src.2013-10-17/src/dplot/create_contour.F
nwchem-6.3.revision2-src.2013-10-17/src/dplot/dplot_dump.F
nwchem-6.3.revision2-src.2013-10-17/src/dplot/dplot.F
nwchem-6.3.revision2-src.2013-10-17/src/dplot/dplot_input.F
nwchem-6.3.revision2-src.2013-10-17/src/dplot/get_transden.F
nwchem-6.3.revision2-src.2013-10-17/src/mcscf/detci/detci_spin.F
nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_analysis.F
nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_davidson.F
nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_init.F
nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_residual.F
nwchem-6.3.revision2-src.2013-10-17/src/nwpw/nwpwlib/Parallel/Parallel-tcgmsg.F
nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_input.F
nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_mo2e_hybrid_2eorb_split.F
nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_mo2e_zones_4a_disk_ga_chop_N5.F
nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_mo2e_zones_4a_disk_ga_N5.F
nwchem-6.3.revision2-src.2013-10-17/src/tools/GNUmakefile
nwchem-6.3.revision2-src.2013-10-17/src/util/util_nwchem_version.F

In other words, there's been changes to the TDDFT module, to the dplot module, TCE  etc.

Looking through the nwchem forum, I think the following posts may hint at what's been changed:

tddft: http://www.nwchem-sw.org/index.php/Special:AWCforum/st/id889/Possible_Bug_in_NWCHEM%3A_TD-B97.html

dplot: http://www.nwchem-sw.org/index.php/Special:AWCforum/st/id1013/Dplot_output_charge_density%2C_tot....html. The integrated (electron) density is printed now.

Not sure about the TCE though.

Here's the (almost fulle) diff -r output:

Only in nwchem-6.3-src.2013-05-28/QA/tests: dplot
diff -r nwchem-6.3-src.2013-05-28/src/config/makefile.h nwchem-6.3.revision2-src.2013-10-17/src/config/makefile.h
2c2
< # $Id: makefile.h 24201 2013-05-09 00:59:44Z edo $
---
> # $Id: makefile.h 24592 2013-09-24 18:49:32Z jhammond $
1171c1171
<         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c22)
---
>         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c22)
1173c1173
<         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c24)
---
>         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c24)
1305c1305
<         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c22)
---
>         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c22)
1307c1307
<         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c24)
---
>         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c24)
1532c1532
<         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c22)
---
>         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c22)
1534c1534
<         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c24)
---
>         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c24)
1697,1700c1697,1704
<        ifdef USE_I4FLAGS
<            ifeq ($(_FC),gfortran)
< #wrong             FOPTIONS += -fdefault-integer-8
<     else  ifeq ($(_FC),crayftn)
---
>       ifeq ($(_FC),gfortran)
>         ifdef USE_I4FLAGS
> #             FOPTIONS += -fdefault-integer-4
>         else
>              FOPTIONS += -fdefault-integer-8
>         endif
>       else ifeq ($(_FC),crayftn)
>         ifdef USE_I4FLAGS
1702,1708c1706
<     else   
<              FOPTIONS += -i4
<            endif
<        else
<          ifeq ($(_FC),gfortran)
<            FOPTIONS += -fdefault-integer-8
<          else  ifeq ($(_FC),crayftn)
---
>         else
1710,1715c1708,1717
<          else
<            FOPTIONS += -i8
<          endif
<        endif
<        DEFINES  += -DEXT_INT
<   MAKEFLAGS = -j 1 --no-print-directory
---
>         endif
>       else
>         ifdef USE_I4FLAGS
>              FOPTIONS += -i4
>         else
>              FOPTIONS += -i8
>         endif
>       endif
>       DEFINES  += -DEXT_INT
>       MAKEFLAGS = -j 1 --no-print-directory
1954c1956
<         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c22)
---
>         GNUMAJOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c22)
1956c1958
<         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null | egrep __VERS | cut -c24)
---
>         GNUMINOR=$(shell $(FC) -dM -E - < /dev/null 2> /dev/null | egrep __VERS | cut -c24)
1969a1972,1974
>         ifeq ($(GNU_GE_4_6),true) 
>           FOPTIMIZE += -march=native -mtune=native
>         else
1974,1976d1978
<         ifeq ($(GNU_GE_4_6),true) 
<           FOPTIMIZE += -march=native -mtune=native
<         else
2198d2199
<    EXPLICITF = TRUE
2211a2213
>     EXPLICITF = TRUE
2224a2227
>     EXPLICITF = TRUE
2247c2250,2251
<     #CC = mpicc
---
>     FC = mpixlf77_r
> 
2248a2253
>         CC         = mpicc
2251,2253c2256,2258
<         FOPTIONS  += -g -funderscoring
<         FOPTIMIZE += -O3 -ffast-math -Wuninitialized 
<         FOPTIMIZE += -O0 -g
---
>         FOPTIONS  += -g -funderscoring -Wuninitialized 
>         FOPTIMIZE += -O3 -ffast-math
>         FDEBUG    += -O1 -g
2262c2267,2281
<         CORE_LIBS +=  -llapack  -lblas 
---
>         CORE_LIBS +=  -llapack $(BLASOPT) -lblas
> 
>         # Here is an example for ALCF:
>         # IBMCMP_ROOT=${IBM_MAIN_DIR}
>         # BLAS_LIB=/soft/libraries/alcf/current/xl/BLAS/lib
>         # LAPACK_LIB=/soft/libraries/alcf/current/xl/LAPACK/lib
>         # ESSL_LIB=/soft/libraries/essl/current/essl/5.1/lib64
>         # XLF_LIB=${IBMCMP_ROOT}/xlf/bg/14.1/bglib64
>         # XLSMP_LIB=${IBMCMP_ROOT}/xlsmp/bg/3.1/bglib64
>         # XLMASS_LIB=${IBMCMP_ROOT}/xlmass/bg/7.3/bglib64
>         # MATH_LIBS="-L${XLMASS_LIB} -lmass -L${LAPACK_LIB} -llapack \
>                      -L${ESSL_LIB} -lesslsmpbg -L${XLF_LIB} -lxlf90_r \
>                      -L${XLSMP_LIB} -lxlsmp -lxlopt -lxlfmath -lxl \
>                      -Wl,--allow-multiple-definition"
>         # Note that ESSL _requires_ USE_64TO32 on Blue Gene
2265,2266d2283
<     FC = mpixlf77_r
<     CC = mpixlc_r
2267a2285,2286
>         EXPLICITF  = TRUE
>         CC         = mpixlc_r
2274,2278d2292
< ifdef USE_I4FLAGS
<         FOPTIONS  = -qintsize=4
< else
<         FOPTIONS  = -qintsize=8 
< endif
2280,2284c2294,2318
<         FOPTIONS  += -qEXTNAME -qxlf77=leadzero
<         FOPTIONS  +=    -qstrict -qthreaded -qnosave -g
<         FOPTIMIZE += -O2 -qarch=qp -qtune=qp -qcache=auto -qunroll=auto -qfloat=rsqrt
< #        FOPTIMIZE += -qhot=level=0 
<         FDEBUG    = -O0 
---
>         ifdef USE_I4FLAGS
>             FOPTIONS = -qintsize=4
>             ifeq ($(BLAS_SIZE),8)
>                 @echo "You cannot use BLAS with 64b integers when"
>                 @echo "the compiler generates 32b integers (USE_I4FLAGS)!"
>                 @exit 1
>             endif # BLAS_SIZE
>         else
>             FOPTIONS = -qintsize=8 
>             ifeq ($(BLAS_SIZE),4)
>                 ifneq ($(USE_64TO32),y)
>                     @echo "You cannot use BLAS with 32b integers when"
>                     @echo "the compiler generates 64b integers unless"
>                     @echo "you do the 64-to-32 conversion!"
>                     @exit 1
>                 endif # USE_64TO32
>             endif # BLAS_SIZE
>         endif # USE_I4FLAGS
> 
>         FDEBUG     = -g -qstrict -O3
>         FOPTIONS  += -g -qEXTNAME -qxlf77=leadzero
>         FOPTIONS  += -qthreaded -qnosave # -qstrict
> #        FOPTIMIZE += -g -O3 -qarch=qp -qtune=qp -qcache=auto -qunroll=auto -qfloat=rsqrt
>         FOPTIMIZE += -O3 -qarch=qp -qtune=qp -qsimd=auto -qhot=level=1 -qprefetch -qunroll=yes #-qnoipa
>         FOPTIMIZE += -qreport -qsource -qlistopt -qlist # verbose compiler output
2425a2460,2466
>   ifeq ($(ARMCI_NETWORK),ARMCI)
>     ifdef EXTERNAL_ARMCI_PATH
>       CORE_LIBS += -L$(EXTERNAL_ARMCI_PATH)/lib -larmci
>     else
>       CORE_LIBS += -L$(NWCHEM_TOP)/src/tools/install/lib -larmci
>     endif
>   else
2426a2468
>   endif
diff -r nwchem-6.3-src.2013-05-28/src/dplot/create_contour.F nwchem-6.3.revision2-src.2013-10-17/src/dplot/create_contour.F
4c4
<      .     no_of_spacings,
---
>      .                          no_of_spacings,tol_rho,
7c7
< * $Id: create_contour.F 19697 2010-10-29 16:57:34Z d3y133 $
---
> * $Id: create_contour.F 24552 2013-08-31 21:23:45Z niri $
30a31
>       double precision tol_rho
46d46
<       Double Precision TOLL
247d246
<          TOLL=1.D-15
257c256
< 
---
> c
259,271c258,271
<      T        TOLL,AO_Bas_Han,g_Dns,
<      &                  nbf_ao_mxnbf_ce,nAtom,1,1,1,
<      U        1,ngrpp,nBF,mBF,.false.,1,
<      &                  Dbl_mb(k_FMat),Dbl_mb(k_PMat),
<      &                  Dbl_mb(k_BMat),0d0,
<      &                  Dbl_mb(k_Scr1),0,0d0,Int_mb(k_ibf),
<      &                  Int_mb(k_iniz),Int_mb(k_ifin),
<      &                  Values(iOffg),0,
<      &              dbl_mb(irchi_atom), 0,
<      &              dbl_mb(k_rdat), int_mb(k_cetobfr),1d0,
<      &               0, 0, .false. )
< 
< 
---
>      T         tol_rho,
>      &         AO_Bas_Han,
>      &         g_Dns,
>      &         nbf_ao_mxnbf_ce,
>      &         nAtom,
>      &         1,1,1,
>      U         1,ngrpp,nBF,mBF,.false.,1,
>      &         Dbl_mb(k_FMat),Dbl_mb(k_PMat),Dbl_mb(k_BMat),0d0,
>      &         Dbl_mb(k_Scr1),0,0d0,Int_mb(k_ibf),
>      &         Int_mb(k_iniz),Int_mb(k_ifin), Values(iOffg),0,
>      &         dbl_mb(irchi_atom),0,
>      &         dbl_mb(k_rdat),int_mb(k_cetobfr),100.d0,
>      &         0, .false. )
> c
diff -r nwchem-6.3-src.2013-05-28/src/dplot/dplot_dump.F nwchem-6.3.revision2-src.2013-10-17/src/dplot/dplot_dump.F
3c3
<      ,     natom,xyz,charge,volume,
---
>      ,     natom,xyz,charge,volume,tol_rho,
17c17
<       double precision spread(3),step(3),angle(3)
---
>       double precision spread(3),step(3),angle(3),tol_rho
34,35c34,35
<             Write(Out_Unit,*)Title
<             Write(Out_Unit,*) 'Total Density'
---
>             Write(Out_Unit,*)"Cube file generated by NWChem"
>             Write(Out_Unit,*) Title
80c80,81
<          if(lgaussian) then
---
> c
>          if(lgaussian) then ! for cube files
85,86c86
<                if(abs(values(i)).lt.1d-10) 
<      .              values(i)=0d0
---
>                if(abs(values(i)).lt.tol_rho) values(i)=0d0
113c113
< c $Id: dplot_dump.F 21176 2011-10-10 06:35:49Z d3y133 $
---
> c $Id: dplot_dump.F 24552 2013-08-31 21:23:45Z niri $
diff -r nwchem-6.3-src.2013-05-28/src/dplot/dplot.F nwchem-6.3.revision2-src.2013-10-17/src/dplot/dplot.F
3c3
< * $Id: dplot.F 24177 2013-05-03 20:42:30Z d3y133 $
---
> * $Id: dplot.F 24552 2013-08-31 21:23:45Z niri $
59a60
>       double precision tol_rho
120a122,125
> c --  Read tol_rho 
>       if (.not. rtdb_get(rtdb, 'dplot:tol_rho', mt_dbl, 1,
>      &   tol_rho)) call errquit('dpinput:rtdbget failed',11, RTDB_ERR)
> c
498a504,505
>           call int_init(rtdb,1,AO_Bas_Han)
>           if (iproc.eq.0) write(luout,*) ' Root: ', iroot
500a508
>           call int_terminate()
563c571
<      .        no_of_spacings,
---
>      .                       no_of_spacings, tol_rho,
607c615
<      ,     natom,dbl_mb(k_xyz),dbl_mb(k_charge),volume,
---
>      ,     natom,dbl_mb(k_xyz),dbl_mb(k_charge),volume,tol_rho,
diff -r nwchem-6.3-src.2013-05-28/src/dplot/dplot_input.F nwchem-6.3.revision2-src.2013-10-17/src/dplot/dplot_input.F
3c3
< * $Id: dplot_input.F 22941 2012-09-30 02:37:23Z niri $
---
> * $Id: dplot_input.F 24552 2013-08-31 21:23:45Z niri $
22c22
<       Parameter (Num_Dirs  = 15)
---
>       Parameter (Num_Dirs  = 16)
40a41
>       double precision tol_rho
45c46
<      A     'dos',
---
>      A     'dos','tol_rho',
61c62,63
<       dodos =.false.
---
>       dodos     =.false.
>       iroot     = 1
103c105
<      &     900, 964, 1997, 9999) ind
---
>      &     900, 964, 1997, 1998, 9999) ind
162d163
<       iroot = 0
252a254
> c
255a258
> c
258a262,269
> c
>  1998 continue
>       tol_rho = 1d-15
>       If (.not. inp_f(tol_rho))
>      &  Call ErrQuit('DPlot_Input: failed to read tol_rho',0,
>      &     INPUT_ERR)
>       goto 10
> c
339a351,356
> *
>       If (.not.rtdb_put(rtdb,'dplot:tol_rho',mt_dbl,
>      &   1,tol_rho))
>      &   Call ErrQuit('DPlot_Input: rtdb_put failed - tol_rho',0,
>      &       RTDB_ERR)
> *
diff -r nwchem-6.3-src.2013-05-28/src/dplot/get_transden.F nwchem-6.3.revision2-src.2013-10-17/src/dplot/get_transden.F
4c4
<      &        g_movecs, g_dens)
---
>      &        g_movecs, g_tdens)
24c24
<          integer basis         ! AO basis set handle
---
>          integer basis            ! AO basis set handle
26c26
<          integer g_dens(ipol)     ! Number of AO basis functions
---
>          integer g_tdens(ipol)    ! Transition density matrix
36c36,37
<          double precision r
---
>          integer icntr,itmom
>          double precision r,cntr(3),tmom(20)
54a56
>          call ga_sync()
58a61
> c        initialization
60c63
<     call ga_zero(g_dens(i))
---
>     call ga_zero(g_tdens(i))
61a65,70
>          do icntr=1,3
>            cntr(icntr)=0.0d0
>          enddo
>          do itmom=1,20
>            tmom(itmom)=0.0d0
>          enddo
77a87,91
>             if (ipol.eq.1) nocc(2)=0
>             if (ipol.eq.1) nmo(2)=0
>             if (ipol.eq.1) nfc(2)=0
>             if (ipol.eq.1) nfv(2)=0
> c
119c133
<            open(unit=69,file=filename,form='formatted',
---
>           open(unit=69,file=filename,form='formatted',
133,135c147,150
<               read(69,*) r  ! energy of root
<               do i=1,ipol
<                if (tda) then
---
>              if (tda) then
>                read(69,*) r  ! energy of root
>                read(69,*) r  ! s2_save(n)
>                do i=1,ipol
140c155,159
<                else
---
>                end do ! ipol
>              else   ! full tddft
>                read(69,*) r  ! energy of root
>                read(69,*) r  ! s2_save(n)
>                do i=1,ipol
144a164,166
>                end do ! ipol
> c
>                do i=1,ipol
149,150c171,172
<                end if  ! tda
<               end do ! ipol
---
>                end do ! ipol
>              end if  ! tda
152c174,175
<            close(unit=69,status='keep',err=1002) ! file
---
>           close(unit=69,status='keep',err=1002) ! file
>           ok = 1
153a177,178
> c
>          call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0)
159c184
<           do i=1,ipol
---
>            do i=1,ipol
162c187
<           enddo
---
>            enddo
169c194,195
<               call ga_copy(g_temp(i),g_dens(i))
---
>           call multipole_density(basis,cntr,3,g_temp(i),tmom,20)  ! transition moments
>           call ga_copy(g_temp(i),g_tdens(i))
174c200,203
<              call tddft_transfm(iroot,g_y,g_movecs,nbf_ao,nocc,nmo,
---
>            do i = 1,ipol
>                 call ga_zero(g_temp(i))
>            end do
>            call tddft_transfm(iroot,g_y,g_movecs,nbf_ao,nocc,nmo,
177,180c206,210
< c            accumulate the Y component of the transition density matrix
<              do i = 1,ipol
<               call ga_add(1.d0,g_dens(i),1.d0,g_temp(i),g_dens(i))
<              end do
---
> c          accumulate the Y component of the transition density matrix
>            do i = 1,ipol
>               call multipole_density(basis,cntr,3,g_temp(i),tmom,20)  ! transition moments
>               call ga_add(1.d0,g_tdens(i),1.d0,g_temp(i),g_tdens(i))
>            end do
182a213,229
>          if (ipol.eq.1) then
>           do i=1,20
>             tmom(i)=tmom(i)*dsqrt(2.0d0)
>           enddo
>          end if 
> c
>          if (ga_nodeid().eq.0) then
>                 write(luout,*) " *** tmom(2)***: ", tmom(2)
>                 write(luout,*) " *** tmom(3)***: ", tmom(3)
>                 write(luout,*) " *** tmom(4)***: ", tmom(4)
>          end if
> c
> c        symmetrize the transition density matrix
>          do i = 1,ipol
>              call ga_symmetrize(g_tdens(i))
>          enddo
> c
186c233
<               Call GA_dAdd(1.d0,g_dens(1),1.d0,g_dens(2),g_dens(1))
---
>               Call GA_dAdd(1.d0,g_tdens(1),1.d0,g_tdens(2),g_tdens(1))
188c235
<               Call GA_dAdd(1.d0,g_dens(1),-1.d0,g_dens(2),g_dens(1))
---
>               Call GA_dAdd(1.d0,g_tdens(1),-1.d0,g_tdens(2),g_tdens(1))
191c238
<                Call GA_Copy(g_dens(2),g_dens(1))
---
>                Call GA_Copy(g_tdens(2),g_tdens(1))
194a242
> c        cleanup
diff -r nwchem-6.3-src.2013-05-28/src/mcscf/detci/detci_spin.F nwchem-6.3.revision2-src.2013-10-17/src/mcscf/detci/detci_spin.F
12c12
< * $Id: detci_spin.F 23708 2013-03-08 21:13:06Z bert $
---
> * $Id: detci_spin.F 24317 2013-06-12 16:58:14Z d3y133 $
162c162,163
<       call ga_access(g_civec, blo, bhi, alo, ahi, k_xxci, bdim)
---
>       if (bhi.gt.0.and.ahi.gt.0) then
>         call ga_access(g_civec, blo, bhi, alo, ahi, k_xxci, bdim)
164c165
< c  Allocate scatter data block and pointer blocks
---
> c       Allocate scatter data block and pointer blocks
166,203c167,206
<       scat_dim = 40000
<       if (.not.ma_push_get(MT_DBL, scat_dim, 'detci:lowdin',
<      $                     l_xa, k_xa))
<      $    call errquit('detci: cannot allocate xa lowdin',0, MA_ERR)
<       if (.not.ma_push_get(MT_INT, scat_dim, 'detci:lowdin',
<      $                     l_ib, k_ib))
<      $    call errquit('detci: cannot allocate ib lowdin',0, MA_ERR)
<       if (.not.ma_push_get(MT_INT, scat_dim, 'detci:lowdin',
<      $                     l_ia, k_ia))
<      $    call errquit('detci: cannot allocate ia lowdin',0, MA_ERR)
<       do istra = alo, ahi
<          call ifill((detci_maxorb*detci_maxorb),0,eij,1)
<          call ifill((detci_maxorb*detci_maxorb),0,pij,1)
<          do iex=1,nexa
<            eij(exa(6,iex,istra),exa(5,iex,istra)) = exa(1,iex,istra)
<            pij(exa(6,iex,istra),exa(5,iex,istra)) = exa(4,iex,istra)
<          enddo
<          offset=(istra-alo)*bdim-blo
<          do istrb = blo, bhi
<             val = -dbl_mb(k_xxci+offset+istrb)
<             if (dabs(val).gt.1.0d-14) then 
<                do iex=1,nexb
<                   iib = exb(5,iex,istrb)
<                   jjb = exb(6,iex,istrb)
<                   if ((eij(iib,jjb).ne.0).and.(pij(iib,jjb).ne.0)) then
<                     jstrb = exb(1,iex,istrb)
<                     jstra = eij(iib,jjb)
<                     xx = val*pij(iib,jjb)*exb(4,iex,istrb)
<                     if (dabs(xx).gt.1.0d-14) then
<                        isc=isc+1
<                        dbl_mb(k_xa+isc-1)=xx
<                        int_mb(k_ib+isc-1)=jstrb
<                        int_mb(k_ia+isc-1)=jstra
<                        if (isc.eq.scat_dim) then
<                           call ga_scatter_acc(g_pvec,dbl_mb(k_xa),
<      &                           int_mb(k_ib),int_mb(k_ia),isc,1.0d0)
<                           isc=0
<                        endif
---
>         scat_dim = 40000
>         if (.not.ma_push_get(MT_DBL, scat_dim, 'detci:lowdin',
>      $                       l_xa, k_xa))
>      $      call errquit('detci: cannot allocate xa lowdin',0, MA_ERR)
>         if (.not.ma_push_get(MT_INT, scat_dim, 'detci:lowdin',
>      $                       l_ib, k_ib))
>      $      call errquit('detci: cannot allocate ib lowdin',0, MA_ERR)
>         if (.not.ma_push_get(MT_INT, scat_dim, 'detci:lowdin',
>      $                       l_ia, k_ia))
>      $      call errquit('detci: cannot allocate ia lowdin',0, MA_ERR)
>         do istra = alo, ahi
>            call ifill((detci_maxorb*detci_maxorb),0,eij,1)
>            call ifill((detci_maxorb*detci_maxorb),0,pij,1)
>            do iex=1,nexa
>              eij(exa(6,iex,istra),exa(5,iex,istra)) = exa(1,iex,istra)
>              pij(exa(6,iex,istra),exa(5,iex,istra)) = exa(4,iex,istra)
>            enddo
>            offset=(istra-alo)*bdim-blo
>            do istrb = blo, bhi
>               val = -dbl_mb(k_xxci+offset+istrb)
>               if (dabs(val).gt.1.0d-14) then 
>                  do iex=1,nexb
>                     iib = exb(5,iex,istrb)
>                     jjb = exb(6,iex,istrb)
>                     if ((eij(iib,jjb).ne.0).and.(pij(iib,jjb).ne.0))
>      &              then
>                       jstrb = exb(1,iex,istrb)
>                       jstra = eij(iib,jjb)
>                       xx = val*pij(iib,jjb)*exb(4,iex,istrb)
>                       if (dabs(xx).gt.1.0d-14) then
>                          isc=isc+1
>                          dbl_mb(k_xa+isc-1)=xx
>                          int_mb(k_ib+isc-1)=jstrb
>                          int_mb(k_ia+isc-1)=jstra
>                          if (isc.eq.scat_dim) then
>                             call ga_scatter_acc(g_pvec,dbl_mb(k_xa),
>      &                             int_mb(k_ib),int_mb(k_ia),isc,1.0d0)
>                             isc=0
>                          endif
>                       endif
205,210c208,212
<                   endif
<                enddo
<             endif
<          enddo
<       enddo
<       if (isc.gt.0) call ga_scatter_acc(g_pvec,dbl_mb(k_xa),
---
>                  enddo
>               endif
>            enddo
>         enddo
>         if (isc.gt.0) call ga_scatter_acc(g_pvec,dbl_mb(k_xa),
212c214,221
<       call ga_release(g_civec, blo, bhi, alo, ahi)
---
>         call ga_release(g_civec, blo, bhi, alo, ahi)
>         if (.not.ma_pop_stack(l_ia))
>      $     call errquit('cannot pop stack ia detci:lowdin',0, MA_ERR)
>         if (.not.ma_pop_stack(l_ib))
>      $     call errquit('cannot pop stack ib detci:lowdin',0, MA_ERR)
>         if (.not.ma_pop_stack(l_xa))
>      $     call errquit('cannot pop stack xa detci:lowdin',0, MA_ERR)
>       endif
214,219d222
<       if (.not.ma_pop_stack(l_ia))
<      $   call errquit('cannot pop stack ia detci:lowdin',0, MA_ERR)
<       if (.not.ma_pop_stack(l_ib))
<      $   call errquit('cannot pop stack ib detci:lowdin',0, MA_ERR)
<       if (.not.ma_pop_stack(l_xa))
<      $   call errquit('cannot pop stack xa detci:lowdin',0, MA_ERR)
diff -r nwchem-6.3-src.2013-05-28/src/nwdft/lr_tddft/tddft_analysis.F nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_analysis.F
6c6
< c $Id: tddft_analysis.F 24091 2013-04-17 17:22:55Z bert $
---
> c $Id: tddft_analysis.F 24553 2013-08-31 21:27:02Z niri $
179a180,182
>       double precision s2_save(nroots) 
>       logical lstores2
>       double precision s2_tmp(nroots)
219c222
< c     CI Vectors file 
---
> c     CI Vectors file
222a226,239
>       if (lcivecs) then
>         do n=1,nroots
>           if (ipol.eq.2) then   ! unrestricted
>             s2_save(n) = 0.0d0
>             s2_tmp(n)  = 0.0d0
>           elseif (singlet) then ! restricted singlets
>             s2_save(n) = 0.0d0
>             s2_tmp(n)  = 0.0d0
>           elseif (triplet) then ! restricted triplets
>             s2_save(n) = 2.0d0
>             s2_tmp(n)  = 2.0d0
>           endif
>         enddo
>       endif
459,494d475
< c --------------------
< c Solution vector file
< c --------------------
< c
<        if (.not.rtdb_cget(rtdb,'tddft:civecs',1,fn_civecs))
<      1  call errquit('tddft_analysis: failed to read vector',0)
< c
<        len_fn_civecs = inp_strlen(fn_civecs)
<        if (singlet) fn_civecs=fn_civecs(1:len_fn_civecs)//"_singlet"
<        if (triplet) fn_civecs=fn_civecs(1:len_fn_civecs)//"_triplet"
< c
<        if (nodezero.and.lcivecs) then
<          write(luout,*) "fn_civecs: ",fn_civecs
<          call util_file_name_resolve(fn_civecs, .false.)
<          open(unit=69,file=fn_civecs,form='formatted',status='unknown')
<          write(LuOut,2010) fn_civecs
<          rewind(69)
<          write(69,*) tda
<          write(69,*) ipol
<          write(69,*) nroots
<          if (ipol.eq.1) nocc(2) = 0
<          write(69,*) nocc(1),nocc(2)
<          if (ipol.eq.1) nmo(2) = 0
<          write(69,*) nmo(1),nmo(2)
<          if (ipol.eq.1) nfc(2) = 0
<          write(69,*) nfc(1),nfc(2)
<          if (ipol.eq.1) nfv(2) = 0
<          write(69,*) nfv(1),nfv(2)
<          if (ipol.eq.1) nov(2) = 0
<          write(69,*) nov(1),nov(2)
<          write(69,*)
<        endif ! nodezero
< c
<  2000 format(/,2x,'No CI vector file is created')
<  2010 format(/,2x,'CI vectors are stored in ',a32)
< c
517,530d497
< c       Write out solution vectors: X (Y=0 in TDA)
< c
<         if (nodezero.and.lcivecs) then
<          do n=1,nroots
<            write(69,*)apbval(n)  ! energy of the root
<            do i=1,ipol
<              do m=1,nov(i)
<                call ga_get(g_x(i),m,m,n,n,r,1)
<                write(69,*) r
<              enddo ! nov
<            enddo ! ipol
<          enddo  ! nroots
<         endif  ! nodezero and lcivecs
< c
557,578d523
< c       g_x = X+Y and g_y = X-Y
< c       Write out vectors: X+Y and X-Y
< c
<         if (nodezero.and.lcivecs) then
<            do n=1,nroots
<              write(69,*)apbval(n) ! energy of the root
<              do i=1,ipol
<                do m=1,nov(i)
<                  call ga_get(g_x(i),m,m,n,n,r,1) ! X vectors
<                  write(69,*) r
<                enddo ! nov
<              enddo ! ipol
< c
<              do i=1,ipol
<                do m=1,nov(i)
<                  call ga_get(g_y(i),m,m,n,n,r,1) ! Y vectors
<                  write(69,*) r
<                enddo ! nov
<              enddo ! ipol
<            enddo ! nroots
<         endif  ! nodezero or lcivecs
< c
588,589d532
<       if (nodezero.and.lcivecs) close(unit=69)
< c
652,653c595
< 
< 
---
> c
877a820
>           if (lcivecs) s2_save(n) = s2
1662a1606,1730
> c
> c ----------------------------------------------------------------------
> c Store the <S2> value for the first cycle of a TDDFT
> c optimization in the RTDB.  This will allow us to use it as a reference
> c for all optimization cycles.
> c ----------------------------------------------------------------------
> c
>       if (lcivecs) then
>         lstores2 = .false.
> c Check if <S2> is already in the RTDB. If it is, we don't do anything
> c else.  Otherwise, we write s2_save to the RTDB.  This only happens if
> c tddft_grad:s2 doesn't exist.
>         if (.not.rtdb_get(rtdb,'tddft_grad:s2',mt_dbl,nroots,s2_tmp))
>      1    lstores2 = .true.
>         if (lstores2) then
>           if (.not.rtdb_put(rtdb,'tddft_grad:s2',mt_dbl,nroots,s2_save))
>      1      call errquit('tddft_analysis: failed to store s2', 0,
>      2        RTDB_ERR)
>         endif
>       endif
> c
> c ---------------------------
> c Handle solution vector file
> c ---------------------------
> c
> c On top of what was present originally for storing
> c the excited state information, we also need <S2> for unrestricted
> c calculations.  This is required because we store every state and
> c it is possible that the states reorder.  We can't use the character
> c of singlet and triplet states to identify states since they can be
> c similar.
> c
>        if (.not.rtdb_cget(rtdb,'tddft:civecs',1,fn_civecs))
>      1  call errquit('tddft_analysis: failed to read vector',0,
>      2    RTDB_ERR)
> c
>        len_fn_civecs = inp_strlen(fn_civecs)
>        if (singlet) fn_civecs=fn_civecs(1:len_fn_civecs)//"_singlet"
>        if (triplet) fn_civecs=fn_civecs(1:len_fn_civecs)//"_triplet"
> c
>        if (nodezero.and.lcivecs) then
>          write(luout,*) "fn_civecs: ",fn_civecs
>          call util_file_name_resolve(fn_civecs, .false.)
>          open(unit=69,file=fn_civecs,form='formatted',status='unknown')
>          write(LuOut,2010) fn_civecs
>          rewind(69)
>          write(69,*) tda
>          write(69,*) ipol
>          write(69,*) nroots
>          if (ipol.eq.1) nocc(2) = 0
>          write(69,*) nocc(1),nocc(2)
>          if (ipol.eq.1) nmo(2) = 0
>          write(69,*) nmo(1),nmo(2)
>          if (ipol.eq.1) nfc(2) = 0
>          write(69,*) nfc(1),nfc(2)
>          if (ipol.eq.1) nfv(2) = 0
>          write(69,*) nfv(1),nfv(2)
>          if (ipol.eq.1) nov(2) = 0
>          write(69,*) nov(1),nov(2)
>          write(69,*)
>        endif ! nodezero
> c
>  2000 format(/,2x,'No CI vector file is created')
>  2010 format(/,2x,'CI vectors are stored in ',a32)
> c
> c ------------
> c Tamm-Dancoff
> c ------------
> c
> c Modified for RPA with B = 0
> c
>       if (tda) then
> c
> c       Write out solution vectors: X (Y=0 in TDA)
> c
>         if (nodezero.and.lcivecs) then
>          do n=1,nroots
>            write(69,*)apbval(n)  ! energy of the root
>            write(69,*)s2_save(n) ! <S2> value of the root
>            do i=1,ipol
>              do m=1,nov(i)
>                call ga_get(g_x(i),m,m,n,n,r,1)
>                write(69,*) r
>              enddo ! nov
>            enddo ! ipol
>          enddo  ! nroots
>         endif  ! nodezero and lcivecs
> c
> c --------------------
> c Full linear response
> c --------------------
> c
>       else  ! full tddft
> c
> c       g_x = X+Y and g_y = X-Y
> c
>         do i=1,ipol
>            call ga_add(1.0d0,g_x(i), 1.0d0,g_y(i),g_x(i)) ! X+Y
>            call ga_add(1.0d0,g_x(i),-2.0d0,g_y(i),g_y(i)) ! X+Y-2Y = X-Y
>         enddo
> c
> c       Write out vectors: X+Y and X-Y
> c
>         if (nodezero.and.lcivecs) then
>            do n=1,nroots
>              write(69,*)apbval(n)  ! energy of the root
>              write(69,*)s2_save(n) ! <S2> value of the root
>              do i=1,ipol
>                do m=1,nov(i)
>                  call ga_get(g_x(i),m,m,n,n,r,1) ! X vectors
>                  write(69,*) r
>                enddo ! nov
>              enddo ! ipol
> c
>              do i=1,ipol
>                do m=1,nov(i)
>                  call ga_get(g_y(i),m,m,n,n,r,1) ! Y vectors
>                  write(69,*) r
>                enddo ! nov
>              enddo ! ipol
>            enddo ! nroots
>         endif  ! nodezero or lcivecs
>       endif ! tda
> c
>       if (nodezero.and.lcivecs) close(unit=69)
diff -r nwchem-6.3-src.2013-05-28/src/nwdft/lr_tddft/tddft_davidson.F nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_davidson.F
8c8
< c $Id: tddft_davidson.F 24076 2013-04-15 16:00:42Z niri $
---
> c $Id: tddft_davidson.F 24309 2013-06-06 18:30:18Z niri $
178d177
<       integer vshift
254,259d252
< c Get reference virtual state
< c --------------------------------------------
<       if (.not.rtdb_get(rtdb,'tddft:vshift',mt_int,1,vshift))
<      &   vshift = 0
< c
< c --------------------------------------------
481c474
<      2          lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2          lowin,owstart,owend,lewin,ewinl,ewinh)
485c478
<      2            lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2            lowin,owstart,owend,lewin,ewinl,ewinh)
500c493
<      2          lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2          lowin,owstart,owend,lewin,ewinl,ewinh)
520c513
<      2            lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2            lowin,owstart,owend,lewin,ewinl,ewinh)
664c657
<      2          lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2          lowin,owstart,owend,lewin,ewinl,ewinh)
668c661
<      2            lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2            lowin,owstart,owend,lewin,ewinl,ewinh)
683c676
<      2          lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2          lowin,owstart,owend,lewin,ewinl,ewinh)
703c696
<      2            lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      2            lowin,owstart,owend,lewin,ewinl,ewinh)
801c794
<      7    diff_max,lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      7    diff_max,lowin,owstart,owend,lewin,ewinl,ewinh)
diff -r nwchem-6.3-src.2013-05-28/src/nwdft/lr_tddft/tddft_init.F nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_init.F
10c10
< c $Id: tddft_init.F 22895 2012-09-23 01:19:55Z niri $
---
> c $Id: tddft_init.F 24357 2013-07-01 22:46:52Z edo $
112a113,114
>       logical xc_got2nd
>       external xc_got2nd
125a128,132
> 
>       if(.not.xc_got2nd()) call errquit(
>      A        'analytic 2nds not ready for these XC functionals',0,
>      &       CAPMIS_ERR)
> 
diff -r nwchem-6.3-src.2013-05-28/src/nwdft/lr_tddft/tddft_residual.F nwchem-6.3.revision2-src.2013-10-17/src/nwdft/lr_tddft/tddft_residual.F
7c7
<      6  diff_max,lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
---
>      6  diff_max,lowin,owstart,owend,lewin,ewinl,ewinh)
9c9
< c $Id: tddft_residual.F 24037 2013-04-11 21:10:58Z bert $
---
> c $Id: tddft_residual.F 24309 2013-06-06 18:30:18Z niri $
88d87
<       integer vshift
828,832c827
<               if (vshift.gt.0) then
<                  k=nocc(i)+1+vshift
<               else
<                  k=mod(l-1,nmo(i)-nfv(i)-nocc(i))+nocc(i)+1
<               end if
---
>               k=mod(l-1,nmo(i)-nfv(i)-nocc(i))+nocc(i)+1
929,933c924
<                 if (vshift.gt.0) then
<                    k=nocc(i)+1+vshift
<                 else
<                    k=mod(l-1,nmo(i)-nfv(i)-nocc(i))+nocc(i)+1
<                 end if
---
>                 k=mod(l-1,nmo(i)-nfv(i)-nocc(i))+nocc(i)+1
diff -r nwchem-6.3-src.2013-05-28/src/nwpw/nwpwlib/Parallel/Parallel-tcgmsg.F nwchem-6.3.revision2-src.2013-10-17/src/nwpw/nwpwlib/Parallel/Parallel-tcgmsg.F
2c2
< * $Id: Parallel-tcgmsg.F 22562 2012-06-05 21:17:04Z bylaska $
---
> * $Id: Parallel-tcgmsg.F 24308 2013-06-06 03:34:42Z jhammond $
894c894
<       /* determine psr - should be made w/o using tmp array! */
---
> c      /* determine psr - should be made w/o using tmp array! */
1033c1033
<       /* determine psr - should be made w/o using tmp array! */
---
> c      /* determine psr - should be made w/o using tmp array! */
diff -r nwchem-6.3-src.2013-05-28/src/tce/tce_input.F nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_input.F
6c6
< c $Id: tce_input.F 24178 2013-05-03 22:05:45Z kowalski $
---
> c $Id: tce_input.F 24360 2013-07-02 18:09:01Z jhammond $
14c14
< c        [FREEZE [[core] (atomic || <integer nfzc default 0>)] \
---
> c        [FREEZE [[core] (atomic || <integer nfzc default 0>)] 
16,18c16,18
< c        [(LCCD||CCD||CCSD||LCCSD||CCSDT||CCSDTQ|| \ 
< c          CCSD(T)||CCSD[T]||QCISD||CISD||CISDT||CISDTQ|| \
< c          MBPT2||MBPT3||MBPT4||MP2||MP3||MP4|| \
---
> c        [(LCCD||CCD||CCSD||LCCSD||CCSDT||CCSDTQ|| 
> c          CCSD(T)||CCSD[T]||QCISD||CISD||CISDT||CISDTQ|| 
> c          MBPT2||MBPT3||MBPT4||MP2||MP3||MP4|| 
diff -r nwchem-6.3-src.2013-05-28/src/tce/tce_mo2e_hybrid_2eorb_split.F nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_mo2e_hybrid_2eorb_split.F
3c3
< C     $Id: tce_mo2e_hybrid_2eorb_split.F 19706 2010-10-29 17:52:31Z d3y133 $
---
> C     $Id: tce_mo2e_hybrid_2eorb_split.F 24292 2013-06-04 01:26:22Z edo $
780c780
<        next = nxtask(-nprocs)
---
>        next = nxtask(-nprocs,1)
1070c1070
<       next = nxtask(-nprocs)
---
>       next = nxtask(-nprocs,1)
1297c1297
<       next = nxtask(-nprocs)
---
>       next = nxtask(-nprocs,1)
1543c1543
<       next = nxtask(-nprocs)
---
>       next = nxtask(-nprocs,1)
1734c1734
<        next = nxtask(-nprocs)
---
>        next = nxtask(-nprocs,1)
diff -r nwchem-6.3-src.2013-05-28/src/tce/tce_mo2e_zones_4a_disk_ga_chop_N5.F nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_mo2e_zones_4a_disk_ga_chop_N5.F
4c4
< C     $Id: tce_mo2e_zones_4a_disk_ga_chop_N5.F 19706 2010-10-29 17:52:31Z d3y133 $
---
> C     $Id: tce_mo2e_zones_4a_disk_ga_chop_N5.F 24330 2013-06-19 22:02:55Z kowalski $
480,483c480,489
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g3b-1),nalength(azone4),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g3b-1),nalength(azone4),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
> c old transpositions
>        call TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & nalength(azone2),nalength(azone1),nalength(azone4),
>      & int_mb(k_range_alpha+g3b-1),
>      & 1,2,4,3,1.0d0)
> c
517,520c523,530
<       CALL TCE_SORT_4KG_(dbl_mb(k_4a),dbl_mb(k_aux),
<      & nalength(azone3),nalength(azone4),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_4a),dbl_mb(k_aux),
> ccx     & nalength(azone3),nalength(azone4),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
>        call TCE_SORT_4(dbl_mb(k_4a),dbl_mb(k_aux),
>      &  nalength(azone2),nalength(azone1),nalength(azone4),
>      &  nalength(azone3),
>      &  1,2,4,3,1.0d0)
553,556c563,572
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g3b-1),nalength(azone3),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g3b-1),nalength(azone3),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
> c  old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      1  nalength(azone2),nalength(azone1),nalength(azone3),
>      2  int_mb(k_range_alpha+g3b-1),
>      3  1,2,4,3,1.0d0)
> c
616,619c632,641
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g3b-1),nalength(azone4),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g3b-1),nalength(azone4),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & nalength(azone2),nalength(azone1),nalength(azone4),
>      & int_mb(k_range_alpha+g3b-1),
>      & 1,2,4,3,1.0d0)
> c
813,816c835,844
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
<      &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
<      &1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
> ccx     &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
> ccx     &1,2,4,3,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & int_mb(k_range_alpha+g2b-1),nalength(azone1),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
> c
853,856c881,889
<       CALL TCE_SORT_4KG_(dbl_mb(k_2g2a),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
<      & nalength(azone1),nalength(azone2), 
<      & 1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_2g2a),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
> ccx     & nalength(azone1),nalength(azone2), 
> ccx     & 1,2,4,3,1.0d0)
> c old transposition
>        CALL TCE_SORT_4(dbl_mb(k_2g2a),dbl_mb(k_aux),
>      & nalength(azone2),nalength(azone1),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
890,893c923,932
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
<      &nalength(azone2),int_mb(k_range_alpha+g2b-1),
<      &1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
> ccx     &nalength(azone2),int_mb(k_range_alpha+g2b-1),
> ccx     &1,2,4,3,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & int_mb(k_range_alpha+g2b-1),nalength(azone2),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
> c
959,962c998,1007
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
<      &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
<      &1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
> ccx     &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
> ccx     &1,2,4,3,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & int_mb(k_range_alpha+g2b-1),nalength(azone1),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
> c
diff -r nwchem-6.3-src.2013-05-28/src/tce/tce_mo2e_zones_4a_disk_ga_N5.F nwchem-6.3.revision2-src.2013-10-17/src/tce/tce_mo2e_zones_4a_disk_ga_N5.F
4c4
< C     $Id: tce_mo2e_zones_4a_disk_ga_N5.F 19706 2010-10-29 17:52:31Z d3y133 $
---
> C     $Id: tce_mo2e_zones_4a_disk_ga_N5.F 24328 2013-06-19 17:52:34Z kowalski $
440,443c440,449
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g3b-1),nalength(azone4),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g3b-1),nalength(azone4),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
> c old transpositions
>        call TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & nalength(azone2),nalength(azone1),nalength(azone4),
>      & int_mb(k_range_alpha+g3b-1),
>      & 1,2,4,3,1.0d0)
> c
477,480c483,491
<       CALL TCE_SORT_4KG_(dbl_mb(k_4a),dbl_mb(k_aux),
<      & nalength(azone3),nalength(azone4),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_4a),dbl_mb(k_aux),
> ccx     & nalength(azone3),nalength(azone4),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
>        call TCE_SORT_4(dbl_mb(k_4a),dbl_mb(k_aux),
>      &  nalength(azone2),nalength(azone1),nalength(azone4),
>      &  nalength(azone3),
>      &  1,2,4,3,1.0d0)
> c
513,516c524,533
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g3b-1),nalength(azone3),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g3b-1),nalength(azone3),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
> c  old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      1  nalength(azone2),nalength(azone1),nalength(azone3),
>      2  int_mb(k_range_alpha+g3b-1),
>      3  1,2,4,3,1.0d0)
> c 
576,579c593,602
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g3b-1),nalength(azone4),
<      & nalength(azone1),nalength(azone2),
<      &2,1,3,4,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g3b-1),nalength(azone4),
> ccx     & nalength(azone1),nalength(azone2),
> ccx     &2,1,3,4,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & nalength(azone2),nalength(azone1),nalength(azone4),
>      & int_mb(k_range_alpha+g3b-1),
>      & 1,2,4,3,1.0d0)
> c
775,778c798,807
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
<      &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
<      &1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
> ccx     &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
> ccx     &1,2,4,3,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & int_mb(k_range_alpha+g2b-1),nalength(azone1),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
> c
815,818c844,852
<       CALL TCE_SORT_4KG_(dbl_mb(k_2g2a),dbl_mb(k_aux),
<      & int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
<      & nalength(azone1),nalength(azone2), 
<      & 1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_2g2a),dbl_mb(k_aux),
> ccx     & int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
> ccx     & nalength(azone1),nalength(azone2), 
> ccx     & 1,2,4,3,1.0d0)
> c old transposition
>        CALL TCE_SORT_4(dbl_mb(k_2g2a),dbl_mb(k_aux),
>      & nalength(azone2),nalength(azone1),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
852,855c886,895
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
<      &nalength(azone2),int_mb(k_range_alpha+g2b-1),
<      &1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
> ccx     &nalength(azone2),int_mb(k_range_alpha+g2b-1),
> ccx     &1,2,4,3,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & int_mb(k_range_alpha+g2b-1),nalength(azone2),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
> c
921,924c961,970
<       CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
<      &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
<      &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
<      &1,2,4,3,1.0d0)
---
> ccx      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_aux),
> ccx     &int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1), 
> ccx     &nalength(azone1),int_mb(k_range_alpha+g2b-1), 
> ccx     &1,2,4,3,1.0d0)
> c old transposition
>       CALL TCE_SORT_4(dbl_mb(k_integral),dbl_mb(k_aux),
>      & int_mb(k_range_alpha+g2b-1),nalength(azone1),
>      & int_mb(k_range_alpha+g3b-1),int_mb(k_range_alpha+g4b-1),
>      & 2,1,3,4,1.0d0)
> c
1011c1057
< c      write(6,*)'DONE --- DONE ---- DONE ---- DONE'
---
> c       write(6,*)'DONE --- DONE ---- DONE ---- DONE'
diff -r nwchem-6.3-src.2013-05-28/src/tools/GNUmakefile nwchem-6.3.revision2-src.2013-10-17/src/tools/GNUmakefile
335a336,338
>     ifdef EXTERNAL_ARMCI_PATH
>         MAYBE_ARMCI = --with-armci=$(EXTERNAL_ARMCI_PATH)
>     else
336a340
>     endif
diff -r nwchem-6.3-src.2013-05-28/src/util/util_nwchem_version.F nwchem-6.3.revision2-src.2013-10-17/src/util/util_nwchem_version.F
4c4
<       nwrev="24277"
---
>       nwrev="24652"
Only in nwchem-6.3-src.2013-05-28/: svnlog

525. Briefly: rotating molecule during optimisation in Gaussian

Most people who have used gaussian will be familiar with molecules that are rotated multiple times during optimisation. Occasionally,  the molecule is rotated repeatedly without any sign of convergence, and without any changes in atom positions.

Uploading the video to blogspot has lead to some distortion (i.e. the molecule is translated -- look at the x,y,z axes in the bottom left) but the general behaviour should still be clear:

The energies aren't really changing:

Solution:
Turning off rotation using 'nosymm' might yield faster convergence:

#P rB3LYP/GEN 5D Opt=()  Freq=() SCF=(MaxCycle=128 ) NoSymm  Punch=(MO)

In ECCE, just uncheck 'Use available symmetry'.

Note that NoSymm doesn't turn off symmetry completely -- it just prevent reorientation:

From http://www.gaussian.com/g_tech/g_ur/k_symmetry.htm

The NoSymmetry keyword prevents molecule reorientation and causes all computations to be performed in the input orientation (although the program still attempts to identify the appropriate point group). Symmetry use can be completely disabled by Symmetry=None; use this option if NoSymm generates an error when identifying the point group.

---

How the figures were made:

First I used this script to extract the .xyz coordinates of the geometry steps:

#!/usr/bin/python
import sys 

def getrawdata(infile):
    f=open(infile,'r')
    opt=0
    geo=0
    struct=[]
    structure=[]
    energies=[]
    energy=[]
    for line in f:
    
        if opt==1 and geo==1 and not ("---" in line):
            structure+=[line.rstrip()]
    
        if 'Coordinates (Angstroms)' in line:
            if opt==0:
                opt=1
                structure=[]
    
        if opt==1 and "--------------------------" in line:
            if geo==0:
                geo=1
            elif geo==1:
                geo=0
                opt=0
        if 'SCF Done' in line:
            energy=filter(None,line.rstrip('\n').split(' ')) 
#       if  'Optimization completed' in line and (opt==0 and geo==0):
            energies+=[float(energy[4])]
            opt=0
            geo=0
            struct+=[structure]
            structure=[]
    
    return struct, energies

def periodictable(elementnumber):
    ptable={1:'H',2:'He',\
    3:'Li', 4:'Be',5:'B',6:'C',7:'N',8:'O',9:'F',10:'Ne',\
    11:'Na',12:'Mg',13:'Al',14:'Si',15:'P',16:'S',17:'Cl',18:'Ar',\
    19:'K',20:'Ca',\
    21:'Sc',22:'Ti',23:'V',24:'Cr',25:'Mn',26:'Fe',27:'Co',28:'Ni',29:'Cu',30:'Zn',\
    31:'Ga',32:'Ge',33:'As',34:'Se',35:'Br',36:'Kr',\
    37:'Rb',38:'Sr',\
    39:'Y',40:'Zr',41:'Nb',42:'Mo',43:'Tc',44:'Ru',45:'Rh',46:'Pd',47:'Ag',48:'Cd',\
    49:'In',50:'Sn',51:'Sb',52:'Te',53:'I',54:'Xe',\
    55:'Cs',56:'Ba',\
    57:'La',58:'Ce',59:'Pr',60:'Nd',61:'Pm',62:'Sm',63:'Eu',64:'Gd',65:'Tb',66:'Dy',67:'Ho',68:'Er',69:'Tm',70:'Yb',71:'Lu',\
    72:'Hf', 73:'Ta', 74:'W',75:'Re', 76:'Os', 77:'Ir',78:'Pt', 79:'Au', 80:'Hg',\
    81:'Tl', 82:'Pb', 83:'Bi',84:'Po',85:'At',86:'Rn',\
    87:'Fr',88:'Ra',\
    89:'Ac',90:'Th',91:'Pa',92:'U',93:'Np',94:'Pu',95:'Am',96:'Cm',97:'Bk',98:'Cf',99:'Es',100:'Fm',101:'Md',102:'No',\
    103:'Lr',104:'Rf',105:'Db',106:'Sg',107:'Bh',108:'Hs',109:'Mt',110:'Ds',111:'Rg',112:'Cn',\
    113:'Uut',114:'Fl',115:'Uup',116:'Lv',117:'Uus',118:'Uuo'}
    element=ptable[elementnumber]
    return element
def genxyzstring(coords,elementnumber):
    x_str='%10.5f'% coords[0]
    y_str='%10.5f'% coords[1]
    z_str='%10.5f'% coords[2]
    element=periodictable(int(elementnumber))
    xyz_string=element+(3-len(element))*' '+10*' '+\
    (8-len(x_str))*' '+x_str+10*' '+(8-len(y_str))*' '+y_str+10*' '+(8-len(z_str))*' '+z_str+'\n'

    return xyz_string

def getstructures(rawdata):

    n=0
    for structure in rawdata:

        n=n+1
        num="%03d" % (n,)
        g=open('structure_'+num+'.xyz','w')
        itson=False
        cartesian=[]

        for item in structure:

            coords=filter(None,item.split(' '))
            coordinates=[float(coords[3]),float(coords[4]),float(coords[5])]
            element=coords[1]
            cartesian+=[genxyzstring(coordinates,element)]
        g.write(str(len(cartesian))+'\n')
        g.write('Structure '+str(n)+'\n')
        for line in cartesian:
            g.write(line)
        g.close()
        cartesian=[]
    return 0

if __name__ == "__main__":
    infile=sys.argv[1]
    rawdata,energies=getrawdata(infile)
    structures=getstructures(rawdata)
    g=open('energies.dat','w')
    for n in range(0,len(energies)):
        g.write(str(n)+'\t'+str(energies[n])+'\n')
    g.close()

Then I turned all the .xyz files into one .xyz file:

cat *.xyz > structures.xyz

Which I then opened in vmd with a view to make a video. Select New Molecule and load structure.xyz. Then go to Extensions, Visualization/Movie Maker. Making a movie using VMD didn't yield to anything of sufficient quality. Instead, I selected Tachyon as the renderer. Under Movie Settings I unchecked Delete Image Files.

Once I had the ppm files I did

ffmpeg -r 2 -i symm.%05d.ppm -vcodec mpeg4 -b 5000k video.avi

The key to getting a decent quality video was picking a high enough bit rate, which is probably obvious to most people who know what bit rate means, but not to a happy amateur like myself (well, now it is).

I made the energy plot using gnuplot and the energies.dat file which the python script above generates.

524. Generating bonds, angles and dihedrals for a molecule for GROMACS

Generating bond, angle and dihedral parameters for GROMACS molecular dynamics simulations is a real PITA when it comes to reasonably large and complex molecules.

Since I am currently looking at the solvation dynamics of a series of isomers of a metal oxide cluster I quickly got tired of working out the bonds and bond constraints for my 101 atom clusters by hand.

So here's a python (2.7.x) script that generates parameters suitable for a GROMACS .top or .itp file.

Let's call the script genparam. You can call it as follows:

genparam molecule.xyz 0.21 3

where molecule.xyz is the molecule in xyz coordinates. 0.21 (nm) is used to generate bonds -- atoms that are 0.21 nm or closer to each other are bonded. The final number should be 1, 2 or 3. If it is 1, then only bonds (and constraints for bonds that are 0.098 nm or shorter -- such as O-H bonds. It's specific for my use, so you can edit the code.) are printed. If it is 2, then bonds and angles are printed. If it is 3, then bonds, angles and dihedrals are printed.

A simple and understandable example using methanol looks like this:

./genparams.py methanol.xyz 0.15 3

[bonds]
[bonds]
1 2 6 0.143 50000.00 ;C OH
1 3 6 0.108 50000.00 ;C H
1 4 6 0.108 50000.00 ;C H
1 5 6 0.108 50000.00 ;C H
[constraints]
2 6 2 0.096 ;OH HO

[ angles ]
2 1 3 1 109.487 50000.00 ;OH C H
2 1 4 1 109.487 50000.00 ;OH C H
2 1 5 1 109.466 50000.00 ;OH C H
1 2 6 1 109.578 50000.00 ;C OH HO
3 1 4 1 109.542 50000.00 ;H C H
3 1 5 1 109.534 50000.00 ;H C H
4 1 5 1 109.534 50000.00 ;H C H
[ dihedrals ]
6 2 1 3 1 60.01 50000.00 2 ;HO OH C H
6 2 1 4 1 60.01 50000.00 2 ;HO OH C H
6 2 1 5 1 180.00 50000.00 2 ;HO OH C H

where methanol.xyz looks like this:

6
Methanol
C      -0.000000010000       0.138569980000       0.355570700000   
OH     -0.000000010000       0.187935770000      -1.074466460000  
H       0.882876920000      -0.383123830000       0.697839450000  
H      -0.882876940000      -0.383123830000       0.697839450000  
H      -0.000000010000       1.145042790000       0.750208830000  
HO     -0.000000010000      -0.705300580000      -1.426986340000

Note that the bond length value may need to be tuned for each molecule -- e.g. 0.21 nm is fine for my metal oxides which have long M-O bonds, while 0.15 nm is better for my methanol.xyz. In the end, you'll probably have to do some manual editing to remove excess bonds, angles and dihedrals.

Also note that there's a switch (prevent_hhbonds) which prevent the formation of bonds between atoms that start with H -- this is to prevent bonds between e.g. H in CH3 units (177 pm). You can change this in the code in the __main__ section.

Finally, note that the function types, multiplicities and forces may need to be edited. I suggest not doing that in the code, but in the output as it will depend on each molecule.

This script won't do science for you, but it'll take care of some of the drudgery of providing geometric descriptors.

Anyway, having said all that, here's the code:

#!/usr/bin/python
import sys
from math import sqrt, pi
from itertools import permutations
from math import acos,atan2

# see table 5.5 (p 132, v 4.6.3) in the gromacs manual
# for the different function types


#from
#http://stackoverflow.com/questions/1984799/cross-product-of-2-different-\
#vectors-in-python
def crossproduct(a, b):
 c = [a[1]*b[2] - a[2]*b[1],
  a[2]*b[0] - a[0]*b[2],
  a[0]*b[1] - a[1]*b[0]]
 return c
#end of code snippet

# mostly from 
# http://www.emoticode.net/python/calculate-angle-between-two-vectors.html
def dotproduct(a,b):
 return sum([a[i]*b[i] for i in range(len(a))])

def veclength(a):
 length=sqrt(a[0]**2+a[1]**2+a[2]**2)
 return length

def ange(a,b,la,lb,angle_unit):
 dp=dotproduct(a,b)
 costheta=dp/(la*lb)
 angle=acos(costheta)
 if angle_unit=='deg':
  angle=angle*180/pi
 return angle
# end of code snippet

def diheder(a,b,c,angle_unit):
 dihedral=atan2(veclength(crossproduct(crossproduct(a,b),\
 crossproduct(b,c))), dotproduct(crossproduct(a,b),crossproduct(b,c)))
 if angle_unit=='deg':
  dihedral=dihedral*180/pi
 return dihedral

def readatoms(infile):
 positions=[]
 f=open(infile,'r')
 atomno=-2
 for line in f:
  atomno+=1
  if atomno >=1:
   position=filter(None,line.rstrip('\n').split(' '))
   if len(position)>3:
    positions+=[[position[0],int(atomno),\
    float(position[1]),float(position[2]),\
    float(position[3])]]
 return positions

def makebonds(positions,rcutoff,prevent_hhbonds):
 
 bonds=[]
 
 for firstatom in positions:
  for secondatom in positions:
   distance=round(sqrt((firstatom[2]-secondatom[2])**2\
+(firstatom[3]-secondatom[3])**2+(firstatom[4]-secondatom[4])**2)/10.0,3)
   xyz=[[firstatom[2],firstatom[3],firstatom[4]],\
[secondatom[2],secondatom[3],secondatom[4]]]
   
   if distance<=rcutoff and firstatom[1]!=secondatom[1]:
    if prevent_hhbonds and (firstatom[0][0:1]=='H' and\
     secondatom[0][0:1]=='H'):
      pass
    elif firstatom[1]<secondatom[1]:
     bonds+=[[firstatom[1],secondatom[1],\
distance,firstatom[0],secondatom[0],xyz[0],xyz[1]]]
    else:
     bonds+=[[secondatom[1],firstatom[1],\
distance,firstatom[0],secondatom[0],xyz[1],xyz[0]]]
 return bonds

def dedupe_bonds(bonds):
 
 newbonds=[]
 for olditem in bonds:
  dupe=False
  for newitem in newbonds:
   if newitem[0]==olditem[0] and newitem[1]==olditem[1]:
    dupe=True
    break;
  if dupe==False:
   newbonds+=[olditem]
 return(newbonds)

def genvec(a,b):
 vec=[b[0]-a[0],b[1]-a[1],b[2]-a[2]]
 return vec

def findangles(bonds,angle_unit):
 # for atoms 1,2,3 we can have the following situations
 # 1-2, 1-3
 # 1-2, 2-3
 # 1-3, 2-3
 # The indices are sorted so that the lower number is always first
 
 angles=[]
 for firstbond in bonds:
    
  for secondbond in bonds:
   if firstbond[0]==secondbond[0] and not (firstbond[1]==\
secondbond[1]): # 1-2, 1-3
    vec=[genvec(firstbond[6],firstbond[5])]
    vec+=[genvec(secondbond[6],secondbond[5])]
    angle=ange(vec[0],vec[1],firstbond[2]*10,\
secondbond[2]*10,angle_unit)
    angles+=[[firstbond[1],firstbond[0],\
secondbond[1],angle,firstbond[4],firstbond[3],secondbond[4],firstbond[6],\
firstbond[5],secondbond[6]]]
   
   if firstbond[0]==secondbond[1] and not (firstbond[1]==\
secondbond[1]): # 1-2, 3-1
    #this should never be relevant since we've sorted the atom numbers
    pass
   
   if firstbond[1]==secondbond[0] and not (firstbond[0]==\
secondbond[1]): # 1-2, 2-3
    vec=[genvec(firstbond[5],firstbond[6])]
    vec+=[genvec(secondbond[6],secondbond[5])]
    angle=ange(vec[0],vec[1],firstbond[2]*10,\
secondbond[2]*10,angle_unit)
    angles+=[[firstbond[0],firstbond[1],\
secondbond[1],angle,firstbond[3],firstbond[4],secondbond[4],firstbond[5],\
firstbond[6],secondbond[6]]]
   
   if firstbond[1]==secondbond[1] and not (firstbond[0]==\
secondbond[0]): # 1-3, 2-3
    vec=[genvec(firstbond[6],firstbond[5])]
    vec+=[genvec(secondbond[6],secondbond[5])]
    angle=ange(vec[0],vec[1],firstbond[2]*10,\
secondbond[2]*10,angle_unit)
    angles+=[[firstbond[0],firstbond[1],\
secondbond[0],angle,firstbond[3],firstbond[4],secondbond[3],firstbond[5],\
firstbond[6],secondbond[5]]]    
 return angles

def dedupe_angles(angles):
 dupe=False
 newangles=[]
 for item in angles:
  dupe=False
  for anotheritem in newangles:
   if item[0]==anotheritem[2] and (item[2]==anotheritem[0]\
 and item[1]==anotheritem[1]):
    dupe=True
    break
  if dupe==False:
   newangles+=[item]
 
 newerangles=[]
 dupe=False
 
 for item in newangles:
  dupe=False
  for anotheritem in newerangles:
   if item[2]==anotheritem[2] and (item[0]==anotheritem[1]\
 and item[1]==anotheritem[0]):
    dupe_item=anotheritem
    dupe=True
    break
  
  if dupe==False:
   newerangles+=[item]
  elif dupe==True:
   if dupe_item[3]>item[3]:
    pass;
   else:
    newerangles[len(newerangles)-1]=item

 newestangles=[]
 dupe=False
 
 for item in newerangles:
  dupe=False
  for anotheritem in newestangles:
   if (sorted(item[0:3]) == sorted (anotheritem[0:3])):
    dupe_item=anotheritem
    dupe=True
    break
  
  if dupe==False:
   newestangles+=[item]
  elif dupe==True:
   if dupe_item[3]>item[3]:
    pass;
   else:
    newestangles[len(newestangles)-1]=item
 return newestangles

def finddihedrals(angles,bonds,angle_unit):
 dihedrals=[]
 for item in angles:
  for anotheritem in bonds:
   if item[2]==anotheritem[0]:
    vec=[genvec(item[7],item[8])]
    vec+=[genvec(item[8],item[9])]
    vec+=[genvec(item[9],anotheritem[6])]
    dihedral=diheder(vec[0],vec[1],vec[2],\
angle_unit)
    dihedrals+=[ [ item[0],item[1],item[2],\
anotheritem[1], dihedral, item[4],item[5],item[6], anotheritem[4] ] ]

   if item[0]==anotheritem[0] and not item[1]==anotheritem[1]:
    vec=[genvec(anotheritem[6],item[7])]
    vec+=[genvec(item[7],item[8])]
    vec+=[genvec(item[8],item[9])]
    dihedral=diheder(vec[0],vec[1],vec[2],\
angle_unit)
    dihedrals+=[ [ anotheritem[1],item[0],item[1],\
item[2], dihedral, anotheritem[4],item[4],item[5], item[6] ] ]
 return dihedrals

def dedup_dihedrals(dihedrals):
 newdihedrals=[]
 
 for item in dihedrals:
  dupe=False
  for anotheritem in newdihedrals:
   rev=anotheritem[0:4]
   rev.reverse()
   if item[0:4] == rev:
    dupe=True
  if not dupe:
   newdihedrals+=[item]
 return newdihedrals

def print_bonds(bonds,func):
 constraints=""
 funcconstr='2'
 print '[bonds]'
 
 for item in bonds:
  if item[2]<=0.098:
   constraints+=(str(item[0])+'\t'+str(item[1])+'\t'+\
funcconstr+'\t'+str("%1.3f"% item[2])+'\t;'+str(item[3])+'\t'+str(item[4])+'\n')
  else: 
   print str(item[0])+'\t'+str(item[1])+'\t'+func+'\t'+\
str("%1.3f"% item[2])+'\t'+'50000.00'+'\t;'+str(item[3])+'\t'+str(item[4])
 print '[constraints]'
 print constraints
 return 0

def print_angles(angles,func):
 print '[ angles ]'
 for item in angles:
  print str(item[0])+'\t'+str(item[1])+'\t'+str(item[2])+'\t'+\
func+'\t'+str("%3.3f" % item[3])+'\t'+'50000.00'+'\t;'\
  +str(item[4])+'\t'+str(item[5])+'\t'+str(item[6])
 return 0

def print_dihedrals(dihedrals,func):
 print "[ dihedrals ]"
 force='50000.00'
 mult='2'
 for item in dihedrals:
  print str(item[0])+'\t'+str(item[1])+'\t'+str(item[2])+'\t'+\
str(item[3])+\
  '\t'+func+'\t'+str("%3.2f" % item[4])+'\t'+force+'\t'+mult+\
  '\t;'+str(item[5])+'\t'+str(item[6])+\
  '\t'+str(item[7])+'\t'+str(item[8])
 return 0

if __name__ == "__main__":
 infile=sys.argv[1]
 rcutoff=float(sys.argv[2]) # in nm
 itemstoprint=int(sys.argv[3])
 
 angle_unit='deg' #{'rad'|'deg'}
 prevent_hhbonds=True # False is safer -- it prevents bonds between
 # atoms whose names start with H

 positions=readatoms(infile)

 bonds=makebonds(positions,rcutoff,prevent_hhbonds)
 bonds=dedupe_bonds(bonds) 
 
 print_bonds(bonds,'6')
 
 angles=findangles(bonds,angle_unit)
 angles=dedupe_angles(angles)
 
 if itemstoprint>=2:
  print_angles(angles,'1')

 dihedrals=finddihedrals(angles,bonds,angle_unit)
 dihedrals=dedup_dihedrals(dihedrals)
 if itemstoprint>=3:
  print_dihedrals(dihedrals,'1')