TABLE OF CONTENTS

::pwtk::calcmol

SYNOPSIS

proc ::pwtk::calcmol {args} {

DESCRIPTION

Utility for making a fast molecule-in-a-box calculation from SMILES specs. Unless explicitly requested otherwise, an orthorhombic "box" will be automatically constructed, based on the size of the molecule.

BEWARE

This routine requires the "obabel" program.

USAGE

    calcmol  ?-relax|-scf?  ?-vacuum VACUUM|-box BOX_SPECS|-nobox?  ?-center|-orig?  \
             ?-mo?  ?-magn NELEC?  ?-smear DEGAUSS?  molecule  SMILES

OPTIONS

Exclusive options: 

    * -vacuum VACUUM ... size of the vacuum in angstroms around the molecule; default = 8 angstroms    
    * -nobox         ... explicitly requests not to assign CELL_PARAMETERS
                         (user must define the lattice parameters prior to calling 'calcmol')    
    * -box BOX_SPECS ... size of the box in angstroms used for CELL_PARAMETERS, where BOX_SPECS are:
       + if BOX_SPECS == A: cubic box (A x A x A),
       + if BOX_SPECS == {A C}: tetragonal box (A x A x C),
       + if BOX_SPECS == {A B C}: orthorhombic (A x B x C).

Non-exclusive options:

ARGUMENTS

Calling ::pwtk::calcmol without options defaults to: -relax -center -vacuum 8.0

RETURN VALUE

The name of the pw.x output file.

SOURCE

    printTitle CALCMOL "Building and calculating a molecular structure from the SMILES specs"
    print CALCMOL "called as:  [procName] $args\n"

    # check if obabel exists
    
    if { [auto_execok obabel] eq "" } {
        ::pwtk::error {"obabel" program does not exist} 1
    }

    # parse options
    
    set options {
        {relax           "make a molecular structural relaxation"}
        {scf             "make a single point SCF calculation"}
        {vacuum.arg 8.0  "size of vacuum around the molecule in angstroms (the same in each direction)"}
        {box.arg    -1   "size of the box in angstroms used for the CELL_PARAMETERS"}
        {nobox           "explicitly requests not to assign CELL_PARAMETERS"}
        {center          "shift the molecule to the center of the simulation BOX"}
        {orig            "shift the molecule to the origin of the simulation BOX"}
        {mo              "calculate signed molecular-orbital densities and write them to XSF file"}
        {magn.arg   0    "total magnetization, i.e., number of unpaired electrons"}
        {smear.arg  -1   "smearing value for occupations = 'smearing' (e.g., use in case of convergency problems)"}
    }
    set usage "?-relax|-scf?  ?-vacuum value | -box {a ?b? ?c?} | -nobox?  ?-center|-orig?  ?-mo?  ?-magn value?  ?-smear value?  molecule  SMILE"

    set narg 2
    parseOpt_
    checkOTypeStrict_ -vacuum $opt(vacuum) double number
    checkOTypeStrict_ -box    $opt(box)    {numberlist double}  "list of numbers"
    checkOTypeStrict_ -magn   $opt(magn)   double number
    checkOTypeStrict_ -smear  $opt(smear)  double number    
    
    # consistency-checks
    
    if { $opt(relax) && $opt(scf) } {
        ::pwtk::error "exclusive -relax and -scf options specified; please specify only one of them" 1
    }
    if { $opt(center) && $opt(orig) } {
        ::pwtk::error "exclusive -center and -orig options specified; please specify only one of them" 1
    }
    
    set lvac [expr {"-vacuum" in $args_}]
    set lbox [expr {"-box" in $args_}]
    if { $lvac + $lbox + $opt(nobox) > 1 } {
        ::pwtk::error "two or more exclusive -vacuum, -box, -nobox options specified; please specify only one of them" 1
    }
    if { $lvac + $lbox + $opt(nobox) == 0 } {
        # -vacuum is the default
        set lvac 1
    }

    set name   [lindex $args 0]
    set smile  [lindex $args 1]
    set prefix [regsub , $name .]; # get rid of ","    
    
    # generate the molecular structure from SMILES
    
    set coor [smilesToCoor $smile]
    set nat  [getNAtoms $coor]
    
    set calc "relax"
    if { $opt(scf) || $nat == 1 } {
        set calc scf
    }
    
    print "Molecule name    :   $name"
    print "SMILES           :   $smile"
    print "Calculation type :   [string toupper $calc]"
    print "Number of atoms  :   $nat\n"
        
    # put a molecule in the box
    
    if { $lvac } {        
        set box [boxify $coor $opt(vacuum)]
        set boxLine "\n   - requested vacuum size (in each direction):  $opt(vacuum) Angstroms"
    } elseif { $lbox } {
        set box [calcmol_boxify_ $opt(box)]
        set boxLine ""
    } else {
        # -nobox
        lassign [calcmol_get_user_box_origin_ $coor] box origin
        set boxLine "\n   - relying on the user defined box"
    }
    set cell [m2cell [matricize [join $box x] ]]

    if { $opt(orig) } {
        set shift origin
        set origin "0.0 0.0 0.0"
    } else {
        ifnotexist origin {
            set origin [lmap o $box {expr 0.5*$o}]
        }
        set shift "center at [eval {format {(%.4f, %.4f, %.4f)}} $origin] Angstroms"
    }
    set coor [eval {::pwtk::atmPosTo $coor} $origin]

    print "MOLECULE in a BOX SETTINGS:${boxLine}
   - box size (in each direction):  [eval {format {%.4f x %.4f x %.4f}} $box] Angstroms
   - molecule shifted to the box $shift\n"

    #
    # make a pw.x calculation
    #
    input_pushpop {
        CONTROL " calculation = [squote $calc], prefix = [squote $prefix] "
        if { $calc == "relax" } {
            IONS ""
        }
        if { ! $opt(nobox) } {
            SYSTEM " ibrav = 0 "
            CELL_PARAMETERS (angstrom) $cell
        }
        ATOMIC_POSITIONS (angstrom) $coor
        ::pwtk::pwi::setNAtoms

        if { [::pwtk::pwi::ibrav] == 0 } {
            print "Box unit-cell vectors:\n[::pwtk::input::cardGet CELL_PARAMETERS]"
        }           
        print "Molecular atomic positions in the box:\n[::pwtk::input::cardGet ATOMIC_POSITIONS]"

        # smearing ?
        if { $opt(smear) > 0.0 } {
            set smearing [::pwtk::input::namelistGetVarValue SYSTEM smearing]
            if { $smearing == {} } { set smearing 'm-p' }
            set occ "occupations = 'smearing', smearing = $smearing, degauss = $opt(smear)"
            print "Smearing requested :"
            print "SYSTEM { $occ }\n"
            SYSTEM $occ
        }

        # set the spin/magnetization setup
        set nbnd [spinSetup $opt(magn) $opt(mo)]

        print CALCMOL "making a pw.x [string toupper $calc] calculation of \"$name\"\n"
        set pwo [runPW pw.$name.$calc]

        # do we want to calculate also MO ?
        
        if { $opt(mo) } {
            printTitle ORBITALS "Calculating valence signed molecular-orbital densities (up to LUMO orbital)"
            INPUTPP "
                filplot = 'psi2.$prefix'
                plot_num  =  7
                kpoint    =  1 
                kband(1)  =  1
                kband(2)  =  $nbnd
                lsign     =  .true.
            "
            PLOT {
                nfile     =  1,
                weight(1) =  1.0,
                iflag     =  3,
                fileout   = '.xsf'
                output_format =  5,
            }
            runPP pp.$name.psi2

            set last [format %03d $nbnd]
            print "Molecular orbital densities written to files :  
   psi2.${prefix}_K001_B001.xsf ... psi2.${prefix}_K001_B$last.xsf\n
To visualize molecular-obital densities, e.g., use:\n
   xcrysden --xsf psi2.${prefix}_K001_B001.xsf\n"
            
            # delete the PP files
            #print "(deleting temporary psi2.${prefix}_K001_B??? files)"
            foreach file [glob -nocomplain psi2.${prefix}_K001_B???] {
                file delete $file
            }
        }
    }
    return $pwo
}