test.py 2.89 KB
Newer Older
Carl Poelking's avatar
Setup.  
Carl Poelking committed
1
2
3
4
5
6
7
8
#! /usr/bin/env python
import soap
import soap.tools

import os
import numpy as np
from momo import osio, endl, flush

9
10
11
12
13
14
15
16
17
18
19
20
21
element_vdw_radius = { 
'C':1.70,
'N':1.55,
'O':1.52,
'H':1.20
}
element_mass = {
'C':12,
'N':14,
'O':16,
'H':1
}

Carl Poelking's avatar
Setup.  
Carl Poelking committed
22
ase_config_list = soap.tools.ase_load_all('configs')
Carl Poelking's avatar
Carl Poelking committed
23
24
for config in ase_config_list:
    print config.config_file
Carl Poelking's avatar
Carl Poelking committed
25
config = ase_config_list[1]
Carl Poelking's avatar
Setup.  
Carl Poelking committed
26
27
osio << config.atoms << endl

Carl Poelking's avatar
Carl Poelking committed
28
29
sigma = 0.5

Carl Poelking's avatar
Setup.  
Carl Poelking committed
30
31
32
# INITIALIZE OPTIONS
options = soap.Options()
options.set('radialbasis.type', 'gaussian')
Carl Poelking's avatar
Carl Poelking committed
33
34
options.set('radialbasis.mode', 'adaptive')
#options.set('radialbasis.mode', 'equispaced')
35
options.set('radialbasis.N', 9)
Carl Poelking's avatar
Carl Poelking committed
36
options.set('radialbasis.sigma', sigma)
37
options.set('radialbasis.integration_steps', 15)
Carl Poelking's avatar
Carl Poelking committed
38
#options.set('radialbasis.N', 9)
Carl Poelking's avatar
Carl Poelking committed
39
options.set('radialcutoff.Rc', 6.8)
Carl Poelking's avatar
Carl Poelking committed
40
41
42
options.set('radialcutoff.Rc_width', 0.5)
options.set('radialcutoff.type', 'shifted-cosine')
options.set('radialcutoff.center_weight', 1.)
Carl Poelking's avatar
Carl Poelking committed
43
#options.set('radialcutoff.center_weight', -7.)
Carl Poelking's avatar
Carl Poelking committed
44

45
options.set('angularbasis.type', 'spherical-harmonic')
46
options.set('angularbasis.L', 6)
Carl Poelking's avatar
Carl Poelking committed
47
48
#options.set('angularbasis.L', 12)
#options.set('densitygrid.N', 20)
49
options.set('densitygrid.N', 20)
Carl Poelking's avatar
Carl Poelking committed
50
options.set('densitygrid.dx', 0.15)
Carl Poelking's avatar
Setup.  
Carl Poelking committed
51
52
53
54

# LOAD STRUCTURE
structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)

55
osio << osio.mg << structure.label << endl
56

57
for atom in structure:
Carl Poelking's avatar
Carl Poelking committed
58
    atom.sigma = sigma # 0.5 # 0.4*element_vdw_radius[atom.type]
59
    atom.weight = 1. #element_mass[atom.type]
Carl Poelking's avatar
Carl Poelking committed
60
61
    #atom.sigma = 0.5*element_vdw_radius[atom.type]
    #atom.weight = element_mass[atom.type]
62
    osio << atom.name << atom.type << atom.weight << atom.sigma << atom.pos << endl
Carl Poelking's avatar
Carl Poelking committed
63
    #if atom.id > 60: atom.weight *= -1
64

Carl Poelking's avatar
Setup.  
Carl Poelking committed
65
66
67
# COMPUTE SPECTRUM
spectrum = soap.Spectrum(structure, options)
spectrum.compute()
68
69
spectrum.computePower()
spectrum.save("test_serialization/%s.spectrum.arch" % structure.label)
Carl Poelking's avatar
Carl Poelking committed
70
spectrum.save("test_invert/%s.spectrum.arch" % structure.label)
71

Carl Poelking's avatar
Carl Poelking committed
72
spectrum.writeDensity(1, "C", "")
Carl Poelking's avatar
Carl Poelking committed
73
74
spectrum.writePowerDensity(1, "C", "", "")
spectrum.writeDensityOnGrid(1, "C", "")
Carl Poelking's avatar
Carl Poelking committed
75
76
77
78
79
#spectrum.writeDensityOnGridInverse(1, "C", "g", "c")

# INVERSION
basis = spectrum.basis
center = spectrum.getAtomic(1, "C").getCenter()
Carl Poelking's avatar
Carl Poelking committed
80
xnkl = spectrum.getAtomic(1, "C").getPower("", "").array
Carl Poelking's avatar
Carl Poelking committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

ofs = open('density.power.python.coeff', 'w')
for n in range(basis.N):
    for k in range(basis.N):
        for l in range(basis.L+1):
            x = xnkl[n*basis.N+k,l]
            ofs.write("%2d %2d %+2d %+1.7e %+1.7e\n" % (n, k, l, x.real, x.imag))
ofs.close()

pox = soap.PowerExpansion(basis)
pox.array = xnkl
ynkl = pox.array

# powex = spectrum.getAtomicSpectrum(1, "C").getPowerExpansion("g", "c")
#osio << osio.mg << "Atomic spectrum" << endl
#atomic = soap.AtomicSpectrum()

Carl Poelking's avatar
Carl Poelking committed
98
99
100
101
#spectrum.writeDensityOnGrid(2, "S", "")
#spectrum.writeDensityOnGrid(7, "C", "") # line.xyz
#spectrum.writeDensityOnGrid(3, "C", "") # linedot.xyz
#spectrum.writeDensityOnGrid(41, "C", "") # C60_pair.xyz
Carl Poelking's avatar
Setup.  
Carl Poelking committed
102
103
104
105
106

osio.okquit()