Inconsistency between high precisions of python and fortran

[Alireza Khorshidi]

This should be closely related to issues #73 and #103, but it is different from issue #103 in the sense that there is no randomness here.

If we run this script:

import numpy as np
from collections import OrderedDict
from ase import Atoms
from ase.calculators.emt import EMT
from amp import Amp
from amp.descriptor.gaussian import Gaussian
from amp.model.neuralnetwork import NeuralNetwork
from amp.model import LossFunction
from amp.regression import Regressor


# The test function for non-periodic systems

convergence = {'energy_rmse': 10.**10.,
               'energy_maxresid': 10.**10.,
               'force_rmse': 10.**10.,
               'force_maxresid': 10.**10., }

regressor = Regressor(optimizer='BFGS')


def non_periodic_0th_bfgs_step_test():

    # Making the list of periodic image

    images = [Atoms(symbols='PdOPd',
                    pbc=np.array([False, False, False], dtype=bool),
                    cell=np.array(
                        [[1.,  0.,  0.],
                         [0.,  1.,  0.],
                            [0.,  0.,  1.]]),
                    positions=np.array(
                        [[0.,  0.,  0.],
                         [0.,  2.,  0.],
                            [0.,  0.,  3.]])),
              Atoms(symbols='Pd2O',
                    pbc=np.array([False, False, False], dtype=bool),
                    cell=np.array(
                        [[1.,  0.,  0.],
                         [0.,  1.,  0.],
                         [0.,  0.,  1.]]),
                    positions=np.array(
                        [[-2., -1., -1.],
                         [1.,  2.,  1.],
                         [3.,  4.,  4.]]))]

    for image in images:
        image.set_calculator(EMT())
        image.get_potential_energy(apply_constraint=False)
        image.get_forces(apply_constraint=False)

    # Parameters

    Gs = {'O': [{'type': 'G2', 'element': 'Pd', 'eta': 0.8},
                {'type': 'G4', 'elements': [
                    'Pd', 'Pd'], 'eta':0.2, 'gamma':0.3, 'zeta':1},
                {'type': 'G4', 'elements': ['O', 'Pd'], 'eta':0.3, 'gamma':0.6,
                 'zeta':0.5}],
          'Pd': [{'type': 'G2', 'element': 'Pd', 'eta': 0.2},
                 {'type': 'G4', 'elements': ['Pd', 'Pd'],
                  'eta':0.9, 'gamma':0.75, 'zeta':1.5},
                 {'type': 'G4', 'elements': ['O', 'Pd'], 'eta':0.4,
                  'gamma':0.3, 'zeta':4}],
                  }

    hiddenlayers = {'O': (2,), 'Pd': (2,),}

    weights = OrderedDict([('O', OrderedDict([(1, np.matrix([[-2.0, 6.0],
                                                             [3.0, -3.0],
                                                             [1.5, -0.90000000000000002],
                                                             [-2.5, -1.5]])),
                                              (2, np.matrix([[5.5],
                                                             [3.6000000000000001,],
                                                             [1.3999999999999999,]]))])),
                           ('Pd', OrderedDict([(1, np.matrix([[-1.0, 3.0],
                                                              [2.0, 4.2000000000000002,],
                                                              [1.0, -0.69999999999999996,],
                                                              [-3.0, 2.0]])),
                                               (2, np.matrix([[4.0],
                                                              [0.5],
                                                              [3.0]]))])),])

    scalings = OrderedDict([('O', OrderedDict([('intercept', -2.2999999999999998,),
                                               ('slope', 4.5)])),
                            ('Pd', OrderedDict([('intercept', 1.6000000000000001,),
                                                ('slope', 2.5)])),])

    for fortran in [True]:
        for cores in range(1, 2):
            label = 'train-nonperiodic/%s-%i' % (fortran, cores)
            print label
            calc = Amp(descriptor=Gaussian(cutoff=6.5,
                                           Gs=Gs,
                                           fortran=fortran,),
                       model=NeuralNetwork(hiddenlayers=hiddenlayers,
                                           weights=weights,
                                           scalings=scalings,
                                           activation='sigmoid',
                                           regressor=regressor,
                                           fortran=fortran,),
                       label=label,
                       dblabel=label,
#                       cores={'localhost': 1},
                        cores=1,
                        )

            lossfunction = LossFunction(convergence=convergence, force_coefficient=0,)
            calc.model.lossfunction = lossfunction
            calc.train(images=images, train_forces=False)

if __name__ == '__main__':
    non_periodic_0th_bfgs_step_test()

with fortran=False, we get the following values for different variables here:

k = 2
float(scaling['slope']) = 2.5
np.dot(np.matrix(ohat[k - 1]).T = [[ 0.047425873177566781058178690955]
 [ 0.425557483188341079127781085845]
 [ 1.                              ]]
np.matrix(delta[k]).T = [[ 0.031179980593771011027071082822]]
dAtomicEnergy_dWeights[symbol][k] = [[ 0.003696844513297943557450508933]
 [ 0.033172185168366272178808173976]
 [ 0.077949951484427529302401183031]]

whereas if run the code with fortran=True, we get the following values for different variables here:

 layer =           2
 p =           1
 q =           1
 unraveled_parameters(element)%slope =   2.5000000000000000     
 ohat(layer)%onedarray(p) =   4.7425873177566781E-002
 delta(layer)%onedarray(q) =   3.1179980593771011E-002
 unraveled_daenergy_dparameters(element)%weights(layer)%twodarray(p, q) =   3.6968445132979431E-003

The last digit of the above value for "unraveled_daenergy_dparameters" is incorrectly manipulated by fortran (!!), and this results in a 10**(-16) difference between fortran and python results for the 20th component of dloss_dparameters vector.