Dihedral calculation

General Discussion
1

Dear all

I did write a matlab program with reads out the dihedral angle from a xyz file as in Avogadro, but in my case I m not able to read out angles bigger then 90 degrees. For example I have an angle (given by avogadro 2: -143) and when i calculate the angle I get from matlab the value of -37 so to get the value from avogadro one has to calculate -180+37, what is the critiria to ajust the value to say when the 180 has to be included in the code, so far I had to look every file and measure it since the wanted value “flips” from around -150 to -62 at some point.

Thank you so much

best jakobc

2

It is not evident if your interest is to access the dihedral / torsion angle as a value, or the generation / improvement of a script / program (which could be implemented in MatLab).

Assuming you aim for the former, and your molecule of interest is encoded as a .sdf file (as below), you can request OpenBabel from the command line to provide a little report about the distances and angles in your model. Enter the folder containing the file in question (perhaps some cd to change into the directory of interest are necessary), then issue a command in pattern of either

obabel -isdf test.sdf -oreport
obabel -isdf test.sdf -oreport > analysis.txt

The first line will provide the results back to the terminal / command line. If you use the second option, the results are redirected (>) into a newly written file, here named analysis.txt in plain ASCII. You can choose freely, the two commands are independent of each other. Below I attach an input file, and the result of the analysis. The latest revision of OpenBabel’s documentation is here.

test.sdf (1.6 KB)
analysis.txt (5.1 KB)

3

Can you share your code and/or the method you’re using?

For example:

Or:

4

Sorry for the delayed answer, here is my matlab code(is is still not working somhow the same issue is still apparing):

function [angle]= Diheadral_angle_reader_one(path,array);

path=[“”];%the path is leading to the file which i will attach at the end of this file, for testing
array=[26 27 29 31]; %orca input starts from 0, the array is incating the number of atom %which is wanted to be concidered in orca numbering! (also for testing)

atom1=array(1)+1;%due to matlab +1, matlab starts to count from one
atom2=array(2)+1;
atom3=array(3)+1;
atom4=array(4)+1;

formatSpeccoord = ‘%s’;% the format
formatSpec = ‘%d’;
fileID = fopen(path);%open xyz
formatSpecatoms=‘%s’;

Sizem = textscan(fileID,formatSpec);
fclose(fileID);
Sizem = uint64(Sizem{1}); %size of molecule that is given in xyz file
C_data0=cell(Sizem,4);
dim=size(C_data0);

fileID = fopen(path);
% C_data0 = textscan(fileID,‘%c %f %f %f’,2,‘Delimiter’,‘\n’);
C = textscan(fileID,formatSpeccoord,‘HeaderLines’,2,‘Delimiter’,‘\n’);
atom1_coord=C{1}{atom1};%char array of the atom type and the respective coorinates
atom2_coord=C{1}{atom2};
atom3_coord=C{1}{atom3};
atom4_coord=C{1}{atom4};

        atom1_coord=str2num(atom1_coord(3:end)); %neglicting the the atom type, 3 because some atom types has 2 letters (eg, vor example Pt) 
        atom2_coord=str2num(atom2_coord(3:end));
        atom3_coord=str2num(atom3_coord(3:end));
        atom4_coord=str2num(atom4_coord(3:end));
        
        vector1=atom1_coord-atom2_coord;% cacluating the vators betweeen the 4 atoms
        vector2=atom2_coord-atom3_coord;
        vector3=atom3_coord-atom4_coord;
        
        
        v1norm=vector1/sqrt(vector1(1)^2+vector1(2)^2+vector1(3)^2); %normalisation of vector1
        v2norm=vector2/sqrt(vector2(1)^2+vector2(2)^2+vector2(3)^2);%normalisation of vector2
        v3norm=vector3/sqrt(vector3(1)^2+vector3(2)^2+vector3(3)^2);%normalisation of vector3
        
        
        normal1=cross(vector1,vector2);% calculating the normal of the vectors
        normal2=cross(vector2,vector3);
        n1=normal1/sqrt(normal1(1)^2+normal1(2)^2+normal1(3)^2);%normal1 but normilsed
        n2=normal2/sqrt(normal1(1)^2+normal1(2)^2+normal1(3)^2);
        
        %i choose the methode of your fist link, with the 16 upvotes,
        %all the other vaiables are named the same as in the link
        m1=cross(normal1,v2norm); %cacluatin the m1 for the first linke, of the fist answer(with 16 upvoes)
        
        x=dot(n1,n2);%
        y=dot(m1,n2);
        control_value=dot(v2norm,n2);% should be 0 and will be a control as said in the post
        
        angle=atan2(y,x);% angle

end

the XYZ file is:

62
Coordinates from ORCA-job 2_opt_S0_DMF_wB97X-D3_def2-TZVP_alpha
C -3.10710226026588 0.24118543720111 -1.11916379438834
C -3.13653697345593 -1.10752175996764 -1.53135088015023
C -4.09370208546276 0.74653577093161 -0.31832217251697
C -4.15316807690156 -1.93334654034725 -1.13854815558684
C -5.18301505643120 -1.43148792088035 -0.31771872237643
C -5.15258052877653 -0.08607603684429 0.09614872978083
N -6.13806605010077 0.41800780498366 0.89344881043909
N -6.20390802745153 -2.25112503383327 0.06268395314644
C -7.12663245937811 -1.75203821303779 0.82248884775704
C -7.08965650301085 -0.38304684387180 1.25528530619125
H -2.34489148028874 -1.48746512050211 -2.16566067103805
H -4.19430217846450 -2.97154307438358 -1.44515385683735
H -2.29252137771917 0.87867529980372 -1.44075820406447
H -4.08744066878187 1.77875014287434 0.01000090485393
S -8.27420198328421 0.23395640220761 2.40521880834832
S -8.37170578750967 -2.85043727524953 1.40540273230530
C -9.79254403646249 -1.80320567852944 1.40735384080210
C -9.75684432458846 -0.53470436593729 1.83323415027339
S -11.16214706727395 0.50600656005298 1.82162244467476
S -11.24025653889001 -2.54664791398125 0.77528214797188
Pt -12.81348919905665 -0.88835944890716 1.05719388305997
S -14.41099296129232 0.77621699209481 1.34357984328353
S -14.47469598372043 -2.24074983916541 0.17518657997524
C -15.80129881453538 0.20568271080354 0.46813824277070
C -15.79703376345837 -1.22710401887076 -0.01956600859505
N -16.85203428590848 0.86117444389771 0.18717531271576
N -16.89961509088588 -1.60211798799019 -0.61189313162456
C -17.13409779848628 -2.86744157653488 -1.29430053580530
H -17.57075431509295 -0.84181657719869 -0.71537615375310
C -17.18033913998781 -4.04186628440231 -0.33474026243415
C -17.76613649384608 -3.93083409260547 0.92048651838374
C -16.68381143005267 -5.27249964855638 -0.74591535567373
C -16.77327940712957 -6.38078703872817 0.08292312748220
C -17.85021157717123 -5.03647121799188 1.75411251555710
C -17.35726537221368 -6.26494534921641 1.33672899143238
H -18.15075533092470 -2.97366413036259 1.25577133941313
H -18.30311816583027 -4.93748636352921 2.73367006005664
H -17.42329324926610 -7.12829102486289 1.98817921039567
H -16.21732285269169 -5.36156508634259 -1.72159135022925
H -16.38163737131382 -7.33469004917215 -0.25026507619352
C -18.43913902975768 -2.75716923730633 -2.07305081142683
H -16.30171623619317 -3.02362348833151 -1.98800281386952
H -19.27599290667424 -2.59796390587548 -1.39014348432072
H -18.61402410904270 -3.68060052581867 -2.62336031498707
H -18.39696532410435 -1.93002784373568 -2.78399018407371
C -16.99268534631188 2.25120697048005 0.57904840094956
C -17.75402702107626 3.00634580895760 -0.49074041013648
C -17.26774414906089 4.21790006991392 -0.96568944639668
C -18.95580593642236 2.51998579090700 -0.99823440942851
C -19.65428138949346 3.23124541026798 -1.96090157653801
C -19.16077055307776 4.44243350850938 -2.43045575180466
C -17.96586142576182 4.93503960314860 -1.92907984649461
H -16.32919682282724 4.60349958172767 -0.58094911537463
H -19.70496032278990 4.99632845354779 -3.18638086769489
H -17.56994604585308 5.87637108157858 -2.29250761259975
H -19.34289346649781 1.57257615481389 -0.64026892539661
H -20.58706876902127 2.83883768226486 -2.34938786704314
C -17.70260825546483 2.33018223684295 1.93219681264278
H -16.00587208076309 2.71645672173098 0.68149837563780
H -18.69302626964000 1.87594750267716 1.86905562778761
H -17.81679416182551 3.37262499457465 2.23288405321831
H -17.12574431098945 1.80768737607124 2.69722819754831

5

With respect, this seems like a bug in your Matlab code. Do you need to do it in Matlab? @Thomas suggested a way to do it using Open Babel’s report, and I posted some Python code.

One thing that’s obvious in your code is that you’re not taking vectors with the correct sign.

vector1 should be atom2_coord - atom1_coord because you’re focusing on the central bond between atom 2 and 3 … so the vector1 needs to point outward from atom 2.

I suspect that’s your problem. Beyond that, I don’t know any matlab, so I’d suggest posting in a forum or finding someone in your group who knows that. (Or switch to Python.)

6

Dear all

I was able to get the dihedral out the correct way (as in avogadro), I did it with the first sugeested link by ghutchris, thank you for giving me the link.

Best
jakobc

7

This topic was automatically closed 3 days after the last reply. New replies are no longer allowed.