remove third_order command and ASE based examples to be added in a new pull request

from ase import Atoms, Atom

from ase.calculators.lammpslib import LAMMPSlib

import numpy as np

import matplotlib.pyplot as plt

from mpi4py import MPI

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

GaAs = Atoms([Atom('Ga', (0.0, 0.0, 0.0)),

              Atom('As', (1.413425, 1.413425, 1.413425))],

             cell=[(0.0, 2.82685, 2.82685), (2.82685, 0.0, 2.82685), (2.82685, 2.82685, 0.0)],

             pbc=True,)

cmds = ["pair_style bop", "pair_coeff * * ../../../../../potentials/GaAs.bop.table Ga As",

        "comm_modify cutoff 12"]

mends = ["info system",

         "dynamical_matrix all eskm 0.000001 file dynmat.dat binary no",

         "neigh_modify delay 0"]

N = 5

GaAs = GaAs.repeat([N, N, N])

lammps = LAMMPSlib(lmpcmds=cmds, atom_types={'Ga': 1, 'As': 2}, amendments=mends, log_file='lammps.log')

GaAs.set_calculator(lammps)

GaAs.get_potential_energy()

if rank == 0:

    dynmat = np.loadtxt("dynmat.dat")

    dynmat = dynmat.reshape(([int(3*(len(dynmat)/3)**0.5), int(3*(len(dynmat)/3)**0.5)]))

    eigv = np.linalg.eigvals(dynmat)

    eigv.sort()

    eigv = np.sqrt(np.abs(eigv))/(2*np.pi)

    plt.hist(eigv, 80)

    plt.xlabel('Frequency (THz)')

    plt.show()

from ase import Atoms, Atom

from ase.calculators.lammpslib import LAMMPSlib

import numpy as np

import matplotlib.pyplot as plt

from mpi4py import MPI

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

GaN = Atoms([Atom('Ga', (1.59, 0.917986928012, 0.0)),

             Atom('Ga', (1.59, -0.917986928012, 2.583)),

             Atom('N', (1.59, 0.917986928012, 1.98891)),

             Atom('N', (1.59, -0.917986928012, 4.57191))],

            cell=[(1.59, -2.75396078403, 0.0), (1.59, 2.75396078403, 0.0), (0.0, 0.0, 5.166)],

            pbc=True)

cmds = ["pair_style tersoff", "pair_coeff * * ../../../../../potentials/GaN.tersoff Ga N"]

mends = ["info system",

         "dynamical_matrix all eskm 0.000001 file dynmat.dat binary no",

         "neigh_modify delay 0"]

N = 6

GaN = GaN.repeat([N, N, N])

lammps = LAMMPSlib(lmpcmds=cmds, atom_types={'Ga': 1, 'N': 2}, amendments=mends, log_file='lammps.log')

GaN.set_calculator(lammps)

GaN.get_potential_energy()

if rank == 0:

    dynmat = np.loadtxt("dynmat.dat")

    dynmat = dynmat.reshape(([int(3*(len(dynmat)/3)**0.5), int(3*(len(dynmat)/3)**0.5)]))

    eigv = np.linalg.eigvals(dynmat)

    eigv.sort()

    eigv = np.sqrt(np.abs(eigv))/(2*np.pi)

    plt.hist(eigv, 80)

    plt.xlabel('Frequency (THz)')

    plt.show()

from ase import Atoms, Atom

from ase.calculators.lammpslib import LAMMPSlib

import numpy as np

import matplotlib.pyplot as plt

from mpi4py import MPI

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

quartz = Atoms(

    [Atom('Si', (1.1545226, -1.99969180169, 0.0)),

     Atom('Si', (1.1545226, 1.99969180169, 3.6036)),

     Atom('Si', (2.6069548, 2.15247249027e-16, 1.8018)),

     Atom('O', (1.6724232, -0.624132037742, 0.64378314)),

     Atom('O', (1.6724232, 0.624132037742, 2.9598186618)),

     Atom('O', (2.1623026, -2.49695388906, 4.2473849418)),

     Atom('O', (3.5392742, 1.13629495821, 1.1580150582)),

     Atom('O', (3.5392742, -1.13629495821, 2.4455813382)),

     Atom('O', (2.1623026, 2.49695388906, 4.76161686))],

    cell=[(2.458, -4.257380885, 0.0), (2.458, 4.257380885, 0.0), (0.0, 0.0, 5.4054)],

    pbc=True,

    )

# number of repeats

N = 3

quartz = quartz.repeat([N, N, N])

header = ['units metal',

          'atom_style charge',

          'atom_modify map array sort 0 0']

cmds = ["pair_style      buck/coul/long 10.0 8.0",

        "pair_coeff	1 1      0      1         0",

        "pair_coeff	1 2   18003.7572 0.20520 133.5381",

        "pair_coeff	2 2    1388.7730 0.36232 175.0000",

        "kspace_style ewald 1.0e-12",

        "set type 1 charge 2.4",

        "set type 2 charge -1.2"]

mends = ["dynamical_matrix all eskm 0.000001 file dynmat.dat binary no",

         "neigh_modify delay 0"]

lammps = LAMMPSlib(lmpcmds=cmds, lammps_header=header, amendments=mends, log_file='lammps.log')

quartz.set_calculator(lammps)

quartz.get_potential_energy()

if rank == 0:

    dynmat = np.loadtxt("dynmat.dat")

    dynmat = dynmat.reshape(([int(3*(len(dynmat)/3)**0.5), int(3*(len(dynmat)/3)**0.5)]))

    eigv = np.linalg.eigvals(dynmat)

    eigv.sort()

    plt.hist(33*np.sqrt(np.abs(eigv))/(2*np.pi), 80)

    plt.xlabel('Frequency (cm-1)')

    plt.show()

# third_order command

## Syntax

```

third_order group-ID style args keyword value ...

```

* group-ID = ID of group of atoms to displace

* style = *regular* or *ballistico*

```

*regular* args = gamma

  gamma = finite difference displacement length

*ballistico* args = gamma

  gamma = finite difference displacement length

```

* zero or more keyword/value pairs may be appended

* keyword = *file* or *binary*

```

*file* value = output_file

  output_file = name of file to dump the dynamical matrix into

*binary* values = *no* or *gzip*

```

## Examples

```

third_order 1 regular 0.000001

third_order 1 ballistico 0.000001

third_order 3 regular 0.00004 file third_order.dat

third_order 5 ballistico 0.00000001 file third_order.dat binary gzip

```

## Description

Calculate the finite difference third order tensor of the selected group.

## Restrictions

None

## Related commands

None

## Default

The option defaults are file = "third_order.dat", binary = no  

# LAMMPS LATTICE DYNAMICS COMMANDS

## THIRD ORDER TENSOR CALCULATOR

This directory contains the ingredients to calculate a third order tensor.  

Example:

```

$THIRD_ORDER=third_order #tensor output file

NP=4 #number of processors

mpirun -np $NP lmp_mpi -in in.silicon -out out.silicon

combine.sh third_order

```

To test out a different silicon example:

```

$THIRD_ORDER=third_order

$LMP_FILE=amorphous_silicon.lmp

cp lmp_bank/$LMP_FILE ./silicon_input_file.lmp

NP=4 #number of processors

mpirun -np $NP lmp_mpi -in in.silicon -out out.silicon

bash combine.sh $THIRD_ORDER

```

## Requires: MANYBODY and MOLECULE packages