test.py 2.9 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
25
for config in ase_config_list:
    print config.config_file
config = ase_config_list[4]
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
73
spectrum.writeDensity(1, "C", "")
spectrum.writePowerDensity(1, "C", "g", "c")
74
#spectrum.writeDensityOnGrid(1, "C", "")
Carl Poelking's avatar
Carl Poelking committed
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#spectrum.writeDensityOnGridInverse(1, "C", "g", "c")

# INVERSION
basis = spectrum.basis
center = spectrum.getAtomic(1, "C").getCenter()
xnkl = spectrum.getAtomic(1, "C").getPower("g", "c").array

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()