def test_ex1():
    import psi4
    import numpy as np
    import resp
    
    # Initialize molecule
    mol = psi4.geometry(""" C   1.45051389  -0.06628932   0.00000000
     H   1.75521613  -0.62865986  -0.87500146
     H   1.75521613  -0.62865986   0.87500146
     H   1.92173244   0.90485897   0.00000000
     C  -0.04233122   0.09849378   0.00000000
     O  -0.67064817  -1.07620915   0.00000000
     H  -1.60837259  -0.91016601   0.00000000
     O  -0.62675864   1.13160510   0.00000000""")
    mol.update_geometry()
    
    # Specify options
    options = {'N_VDW_LAYERS'       : 4,
               'VDW_SCALE_FACTOR'   : 1.4,
               'VDW_INCREMENT'      : 0.2,
               'VDW_POINT_DENSITY'  : 1.0,
               'resp_a'             : 0.0005,
               'RESP_B'             : 0.1,
               }
    
    # Call for first stage fit
    charges1 = resp.resp([mol], [options])
    print('Electrostatic Potential Charges')
    print(charges1[0][0])
    print('Restrained Electrostatic Potential Charges')
    print(charges1[0][1])
    # Reference charges are generated by the R.E.D.-III.5 tools
    # with GAMESS as the quantum chemistry package
    reference_charges1 = np.array([-0.294974,  0.107114,  0.107114,  0.084795,
                                    0.803999, -0.661279,  0.453270, -0.600039])
    print('Reference RESP Charges')
    print(reference_charges1)
    print('Difference')
    print(charges1[0][1]-reference_charges1)

    assert np.allclose(charges1[0][1], reference_charges1, atol=5e-4)
    
    # Change the value of the RESP parameter A
    options['resp_a'] = 0.001
    
    # Add constraint for atoms fixed in second stage fit
    constraint_charge = []
    for i in range(4, 8):
        constraint_charge.append([charges1[0][1][i], [i+1]])
    options['constraint_charge'] = constraint_charge
    options['constraint_group'] = [[2, 3, 4]]
    options['grid'] = '1_%s_grid.dat' %mol.name()
    options['esp'] = '1_%s_grid_esp.dat' %mol.name()
    mol.set_name('stage2')
    
    # Call for second stage fit
    charges2 = resp.resp([mol], [options])
    
    # Get RESP charges
    print("\nStage Two:\n")
    print('RESP Charges')
    print(charges2[0][1])
    reference_charges2 = np.array([-0.290893,  0.098314,  0.098314,  0.098314,
                                   0.803999, -0.661279,  0.453270, -0.600039])
    print('Reference RESP Charges')
    print(reference_charges2)
    print('Difference')
    print(charges2[0][1]-reference_charges2)

    assert np.allclose(charges2[0][1], reference_charges2, atol=5e-4)
