fixed a bug in building the dihedral list

[zwang123]

Hi,

In the previous one, it is possible to find dihedrals with repeating atoms, leading to DivisionByZero error (i.e. generating nan). The fixed code resolved this issue.

I do not know if we need to modify angles_r = [Angle(center[-2], center[-1], j).eval(coords) for j in neigh_r] to ... center[0], center[-1] ... Please modify it accordingly if necessary.

[jhrmnn]

Thanks! And sorry for a late response. Could you please give an example of a geometry that was failing and that this patch fixes? I'll include it in the test suite.

[zwang123]

Please try this one. It was done a while ago, so I am not sure if this was my geometry that failed using the original codes.

O 7.098637823416 16.971023444792 6.380247074350 H 7.197946705598 17.877545936880 6.714905976001 H 7.914514819435 16.534616934031 5.908851403895 H 6.088748006295 16.403365566544 6.560590428393 O 4.595372130861 14.424819756679 9.021009962591 H 5.236106969531 14.003807451813 9.640189273152 H 3.678351107573 14.390572637584 9.375270964551 O 5.020575188210 15.791160980651 6.758521185763 H 4.314164242355 15.763235277936 6.085426399656 H 4.857492840207 15.273171434313 7.612775348354 O 9.055933636453 15.904286482443 5.244078633375 H 9.270306093544 16.040547601426 4.309397172129 H 9.736390436523 15.367147494908 5.690336177789

[jhrmnn]

I've tried with a version from 12 Jul, and I'm getting no atom index repetitions in the dihedrals for this geometry.

coords = InternalCoords(geom)
assert not [dih for dih in coords.dihedrals if len(set(dih.idx)) < 4]

Could you try to find the geometry for which this was happening? I was trying to understand how could the situation with the repetitions occur, and I just cannot imagine it.

[zwang123]

Maybe try this one: New geometry (unit Bohr) 1 H -0.000000000000 0.000000000000 -1.142569988888 2 O 1.784105551801 1.364934064507 -1.021376180623 3 H 2.248320553963 2.318104360291 -2.500037742933 4 H 3.285761299420 0.674554743661 -0.259576564237 5 O -1.784105551799 -1.364934064536 -1.021376180591 6 H -2.248320553963 -2.318104360291 -2.500037742933 7 H -3.285761299424 -0.674554743614 -0.259576564287 8 O 5.839754502206 -0.500682935209 1.037064691223 9 H 7.440059622286 -1.597667062287 0.565115038647 10 H 6.475526400773 0.638572472561 2.500357106648 11 O -5.839754502205 0.500682935191 1.037064691242 12 H -7.440059622286 1.597667062287 0.565115038647 13 H -6.475526400773 -0.638572472561 2.500357106648

[zwang123]

Below is my code with pyscf 1.5.3+pyberny 0.3.2:

from pyscf import gto, dft, scf, cc, geomopt, lib

mol = gto.M(
          atom=read_coords(xyzname, format="xyz") # This is my function that reads xyz file
        , charge=1
        , verbose=4
        , basis='6311g')

mf = dft.RKS(mol)
mf.xc = 'BLYP'
opt_mol = geomopt.optimize(mf, assert_convergence=True, verbose=4)

[jhrmnn]

Thanks! I can reproduce it now. And I understand the first part of the fix, I didn't think about weak dihedrals.

For the second part, I just realized that the condition should be perhaps much more stringent? Already now, only angles larger than 45 degrees are added, so perhaps also only dihedrals where the angles are larger than 45 degrees should be added. So the check should really be ang >= pi/4 rather than just ang >= 1e-3.

Anyway, thanks for the fix.