TABLE OF CONTENTS

ABINIT/calc_b_matrix [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_b_matrix

FUNCTION

  calculate values of derivatives of internal coordinates as a function of
  cartesian ones =  B matrix

COPYRIGHT

 Copyright (C) 2003-2018 ABINIT group (MVer)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

 angs= number of angles
 bonds(2,2,nbond)=for a bond between iatom and jatom
              bonds(1,1,nbond) = iatom
              bonds(2,1,nbond) = icenter
              bonds(1,2,nbond) = jatom
              bonds(2,2,nbond) = irshift
 carts(2,ncart)= index of total primitive internal, and atom (carts(2,:))
 dihedrals(2,4,ndihed)=indexes to characterize dihedrals
 nang(2,3,nang)=indexes to characterize angles
 nbond=number of bonds
 ncart=number of auxiliary cartesian atom coordinates (used for constraints)
 ndihed= number of dihedrals
 ninternal=nbond+nang+ndihed+ncart: number of internal coordinates
 nrshift= dimension of rshift
 rprimd(3,3)=dimensional real space primitive translations (bohr)
 rshift(3,nrshift)=shift in xred that must be done to find all neighbors of
                   a given atom within a given number of neighboring shells
 xcart(3,natom)=cartesian coordinates of atoms (bohr)

OUTPUT

 b_matrix(ninternal,3*natom)=matrix of derivatives of internal coordinates
   wrt cartesians

SIDE EFFECTS

NOTES

PARENTS

      xcart2deloc

CHILDREN

      acrossb

SOURCE

 53 #if defined HAVE_CONFIG_H
 54 #include "config.h"
 55 #endif
 56 
 57 subroutine calc_b_matrix(deloc,natom,rprimd,xcart,b_matrix)
 58 
 59  use defs_basis
 60  use m_abimover
 61  use m_profiling_abi
 62 
 63 !This section has been created automatically by the script Abilint (TD).
 64 !Do not modify the following lines by hand.
 65 #undef ABI_FUNC
 66 #define ABI_FUNC 'calc_b_matrix'
 67  use interfaces_45_geomoptim, except_this_one => calc_b_matrix
 68 !End of the abilint section
 69 
 70  implicit none
 71 
 72 !Arguments ------------------------------------
 73 !scalars
 74  integer,intent(in) :: natom
 75  type(delocint),intent(inout) :: deloc
 76 
 77 !arrays
 78  real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
 79  real(dp),intent(out) :: b_matrix(deloc%ninternal,3*natom)
 80 
 81 !Local variables-------------------------------
 82 !scalars
 83  integer :: i1,i2,i3,i4,iang,ibond,icart,idihed,iprim,s1,s2,s3,s4
 84 !arrays
 85  real(dp) :: bb(3),r1(3),r2(3),r3(3),r4(3)
 86 
 87 ! *************************************************************************
 88 
 89  iprim=0
 90  b_matrix(:,:) = zero
 91 
 92  do ibond=1,deloc%nbond
 93    i1 = deloc%bonds(1,1,ibond)
 94    s1 = deloc%bonds(2,1,ibond)
 95    r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
 96 &   +deloc%rshift(2,s1)*rprimd(:,2)&
 97 &   +deloc%rshift(3,s1)*rprimd(:,3)
 98    i2 = deloc%bonds(1,2,ibond)
 99    s2 = deloc%bonds(2,2,ibond)
100    r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
101 &   +deloc%rshift(2,s2)*rprimd(:,2)&
102 &   +deloc%rshift(3,s2)*rprimd(:,3)
103    iprim=iprim+1
104    call dbond_length_d1(r1,r2,bb)
105    b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
106    call dbond_length_d1(r2,r1,bb)
107    b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
108  end do
109 
110 !second: angle values (ang)
111  do iang=1,deloc%nang
112    i1 = deloc%angs(1,1,iang)
113    s1 = deloc%angs(2,1,iang)
114    r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
115 &   +deloc%rshift(2,s1)*rprimd(:,2)&
116 &   +deloc%rshift(3,s1)*rprimd(:,3)
117    i2 = deloc%angs(1,2,iang)
118    s2 = deloc%angs(2,2,iang)
119    r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
120 &   +deloc%rshift(2,s2)*rprimd(:,2)&
121 &   +deloc%rshift(3,s2)*rprimd(:,3)
122    i3 = deloc%angs(1,3,iang)
123    s3 = deloc%angs(2,3,iang)
124    r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
125 &   +deloc%rshift(2,s3)*rprimd(:,2)&
126 &   +deloc%rshift(3,s3)*rprimd(:,3)
127    iprim=iprim+1
128    call dang_d1(r1,r2,r3,bb)
129    b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
130    call dang_d2(r1,r2,r3,bb)
131    b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
132    call dang_d1(r3,r2,r1,bb)
133    b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:)
134  end do
135 
136 !third: dihedral values
137  do idihed=1,deloc%ndihed
138    i1 = deloc%dihedrals(1,1,idihed)
139    s1 = deloc%dihedrals(2,1,idihed)
140    r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
141 &   +deloc%rshift(2,s1)*rprimd(:,2)&
142 &   +deloc%rshift(3,s1)*rprimd(:,3)
143    i2 = deloc%dihedrals(1,2,idihed)
144    s2 = deloc%dihedrals(2,2,idihed)
145    r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
146 &   +deloc%rshift(2,s2)*rprimd(:,2)&
147 &   +deloc%rshift(3,s2)*rprimd(:,3)
148    i3 = deloc%dihedrals(1,3,idihed)
149    s3 = deloc%dihedrals(2,3,idihed)
150    r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
151 &   +deloc%rshift(2,s3)*rprimd(:,2)&
152 &   +deloc%rshift(3,s3)*rprimd(:,3)
153    i4 = deloc%dihedrals(1,4,idihed)
154    s4 = deloc%dihedrals(2,4,idihed)
155    r4(:) = xcart(:,i4)+deloc%rshift(1,s4)*rprimd(:,1)&
156 &   +deloc%rshift(2,s4)*rprimd(:,2)&
157 &   +deloc%rshift(3,s4)*rprimd(:,3)
158 !  write(std_out,*) 'dihed ',idihed
159 !  write(std_out,*) r1
160 !  write(std_out,*) r2
161 !  write(std_out,*) r3
162 !  write(std_out,*) r4
163 
164    iprim=iprim+1
165    call ddihedral_d1(r1,r2,r3,r4,bb)
166    b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
167    call ddihedral_d2(r1,r2,r3,r4,bb)
168    b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
169    call ddihedral_d2(r4,r3,r2,r1,bb)
170    b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:)
171    call ddihedral_d1(r4,r3,r2,r1,bb)
172    b_matrix(iprim,3*(i4-1)+1:3*i4) = b_matrix(iprim,3*(i4-1)+1:3*i4) + bb(:)
173  end do
174 
175  do icart=1,deloc%ncart
176    iprim=iprim+1
177    b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) = &
178 &   b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) + one
179  end do
180 
181 !DEBUG
182 ! write (200,*) 'calc_b_matrix : b_matrix = '
183 ! do iprim=1,deloc%ninternal
184 !   do i1=1, 3*natom
185 !     write (200,'(E16.6,2x)',ADVANCE='NO') b_matrix(iprim,i1)
186 !   end do
187 !   write (200,*)
188 ! end do
189 !ENDDEBUG
190 
191 
192 end subroutine calc_b_matrix

ABINIT/dang_d1 [ Functions ]

[ Top ] [ Functions ]

NAME

 dang_d1

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

257 subroutine dang_d1(r1,r2,r3,bb)
258 
259  use defs_basis
260  use m_profiling_abi
261  use m_abimover
262 
263  use m_geometry,  only : acrossb
264 
265 !This section has been created automatically by the script Abilint (TD).
266 !Do not modify the following lines by hand.
267 #undef ABI_FUNC
268 #define ABI_FUNC 'dang_d1'
269 !End of the abilint section
270 
271  implicit none
272 
273 !Arguments ------------------------------------
274 !arrays
275  real(dp),intent(in) :: r1(3),r2(3),r3(3)
276  real(dp),intent(out) :: bb(3)
277 
278 !Local variables ------------------------------
279 !scalars
280  real(dp) :: cos_ang,n1,n1232,n2,tmp
281 !arrays
282  real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3)
283 
284 !************************************************************************
285  n1=bond_length(r1,r2)
286  n2=bond_length(r3,r2)
287 
288  rpt12(:) = r1(:)-r2(:)
289  rpt32(:) = r3(:)-r2(:)
290 
291  cos_ang = (rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3))/n1/n2
292  if (cos_ang > one - epsilon(one)*two) then
293    cos_ang = one
294  else if(cos_ang < -one + epsilon(one)*two) then
295    cos_ang = -one
296  end if
297 
298  rpt(:) = rpt32(:)/n1/n2 - rpt12(:)*cos_ang/n1/n1
299 
300  tmp = sqrt(one-cos_ang**2)
301  bb(:) = zero
302  if (tmp > epsilon(one)) then
303    bb(:) = rpt(:) * (-one)/tmp
304  end if
305 
306 !TEST: version from MOLECULAR VIBRATIONS EB Wilson
307  call acrossb(rpt12,rpt32,cp1232)
308  n1232 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
309  rpt(:) = (cos_ang*rpt12(:)*n2/n1 - rpt32(:))/n1232
310  if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then
311    write(std_out,*) 'Compare bb ang 1 : '
312    write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
313  end if
314  bb(:) = rpt(:)
315 
316 end subroutine dang_d1

ABINIT/dang_d2 [ Functions ]

[ Top ] [ Functions ]

NAME

 dang_d2

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

335 subroutine dang_d2(r1,r2,r3,bb)
336 
337  use defs_basis
338  use m_profiling_abi
339  use m_abimover
340 
341  use m_geometry,  only : acrossb
342 
343 !This section has been created automatically by the script Abilint (TD).
344 !Do not modify the following lines by hand.
345 #undef ABI_FUNC
346 #define ABI_FUNC 'dang_d2'
347 !End of the abilint section
348 
349  implicit none
350 
351 !Arguments ------------------------------------
352 !arrays
353  real(dp),intent(in) :: r1(3),r2(3),r3(3)
354  real(dp),intent(out) :: bb(3)
355 
356 !Local variables ------------------------------
357 !scalars
358  real(dp) :: cos_ang,n1,n1232,n2,tmp
359 !arrays
360  real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3)
361 
362 !************************************************************************
363  n1=bond_length(r1,r2)
364  n2=bond_length(r3,r2)
365 
366  rpt12(:) = r1(:)-r2(:)
367  rpt32(:) = r3(:)-r2(:)
368 
369  cos_ang = (rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3))/n1/n2
370  if (cos_ang > one - epsilon(one)*two) then
371    cos_ang = one
372  else if(cos_ang < -one + epsilon(one)*two) then
373    cos_ang = -one
374  end if
375 
376  rpt(:) = -rpt32(:)/n1/n2 - rpt12(:)/n1/n2 &
377 & + rpt12(:)*cos_ang/n1/n1 + rpt32(:)*cos_ang/n2/n2
378 
379  tmp = sqrt(one-cos_ang**2)
380  bb(:) = zero
381  if (tmp > tol12) then
382    bb(:) = rpt(:) * (-one)/tmp
383  end if
384 
385 !TEST: version from MOLECULAR VIBRATIONS EB Wilson
386  call acrossb(rpt12,rpt32,cp1232)
387  n1232 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
388  rpt(:) = ((n1-n2*cos_ang)*rpt12(:)/n1 + (n2-n1*cos_ang)*rpt32(:)/n2) / n1232
389  if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3))  > tol10) then
390    write(std_out,*) 'Compare bb ang 2 : '
391    write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
392  end if
393  bb(:) = rpt(:)
394 
395 
396 end subroutine dang_d2

ABINIT/dbond_length_d1 [ Functions ]

[ Top ] [ Functions ]

NAME

 dbond_length_d1

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

211 subroutine dbond_length_d1(r1,r2,bb)
212 
213  use defs_basis
214  use m_profiling_abi
215  use m_abimover
216 
217 !This section has been created automatically by the script Abilint (TD).
218 !Do not modify the following lines by hand.
219 #undef ABI_FUNC
220 #define ABI_FUNC 'dbond_length_d1'
221 !End of the abilint section
222 
223  implicit none
224 
225 !Arguments ------------------------------------
226 !arrays
227  real(dp),intent(in) :: r1(3),r2(3)
228  real(dp),intent(out) :: bb(3)
229 
230 !Local variables ------------------------------
231 !arrays
232  real(dp) :: rpt(3)
233 
234 !************************************************************************
235  rpt(:) = r1(:)-r2(:)
236  bb(:) = rpt(:)/bond_length(r1,r2)
237 
238 end subroutine dbond_length_d1

ABINIT/ddihedral_d1 [ Functions ]

[ Top ] [ Functions ]

NAME

 ddihedral_d1

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

415 subroutine ddihedral_d1(r1,r2,r3,r4,bb)
416 
417  use defs_basis
418  use m_profiling_abi
419  use m_geometry,  only : acrossb
420 
421 !This section has been created automatically by the script Abilint (TD).
422 !Do not modify the following lines by hand.
423 #undef ABI_FUNC
424 #define ABI_FUNC 'ddihedral_d1'
425 !End of the abilint section
426 
427  implicit none
428 
429 !Arguments ------------------------------------
430 !arrays
431  real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3)
432  real(dp),intent(out) :: bb(3)
433 
434 !Local variables ------------------------------------
435 !scalars
436  real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,tmp
437 !arrays
438  real(dp) :: cp1232(3),cp32_1232(3),cp32_3432(3),cp3432(3),cpcp(3),rpt(3)
439  real(dp) :: rpt12(3),rpt32(3),rpt34(3)
440 
441 !******************************************************************
442  rpt12(:) = r1(:)-r2(:)
443  rpt32(:) = r3(:)-r2(:)
444  rpt34(:) = r3(:)-r4(:)
445 
446  call acrossb(rpt12,rpt32,cp1232)
447  call acrossb(rpt34,rpt32,cp3432)
448 
449 !DEBUG
450 !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232
451 !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432
452 !ENDDEBUG
453 
454  n1 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
455  n2 = sqrt(cp3432(1)**2+cp3432(2)**2+cp3432(3)**2)
456 
457  cos_dihedral = (cp1232(1)*cp3432(1)+cp1232(2)*cp3432(2)+cp1232(3)*cp3432(3))/n1/n2
458  if (cos_dihedral > one - epsilon(one)*two) then
459    cos_dihedral = one
460  else if(cos_dihedral < -one + epsilon(one)*two) then
461    cos_dihedral = -one
462  end if
463 !we use complementary of standard angle, so
464 !cos_dihedral = -cos_dihedral
465 
466  call acrossb(cp1232,cp3432,cpcp)
467  cpcp(:) = cpcp(:)/n1/n2
468 !we use complementary of standard angle, but sin is invariant
469  sin_dihedral = -(cpcp(1)*rpt32(1)+cpcp(2)*rpt32(2)+cpcp(3)*rpt32(3))&
470 & /sqrt(rpt32(1)**2+rpt32(2)**2+rpt32(3)**2)
471  dih_sign = one
472  if (sin_dihedral < -epsilon(one)) then
473    dih_sign = -one
474  end if
475 
476 !DEBUG
477 !write(std_out,'(a,3E16.6)') 'ddihedral_d1 : cos abs(sin) dih_sign= ',&
478 !&    cos_dihedral,sin_dihedral,dih_sign
479 !ENDDEBUG
480 
481 !ddihedral_d1 = dih_sign* acos(cos_dihedral)
482  call acrossb(rpt32,cp1232,cp32_1232)
483  call acrossb(rpt32,cp3432,cp32_3432)
484 
485  rpt(:) = cp32_3432(:)/n1/n2 - cp32_1232(:)/n1/n1 * cos_dihedral
486  bb(:) = zero
487 
488 !DEBUG
489 !write(std_out,*) 'ddihedral_d1 cp1232 cp3432 = ',cp1232,cp3432,rpt32
490 !write(std_out,*) 'ddihedral_d1 cp32_1232 cp32_3432 = ',cp32_1232,cp32_3432,cos_dihedral,n1,n2
491 !write(std_out,*) 'ddihedral_d1 rpt = ',rpt
492 !ENDDEBUG
493 
494  tmp = sqrt(one-cos_dihedral**2)
495  if (tmp > tol12) then
496 !  we use complementary of standard angle, so cosine in acos has - sign,
497 !  and it appears for the derivative
498    bb(:) = -dih_sign * rpt(:) * (-one) / tmp
499  else
500    bb(:) = dih_sign * cp32_3432(:) / n1 / n2 / &
501 &   sqrt(cp32_3432(1)**2+cp32_3432(2)**2+cp32_3432(3)**2)
502  end if
503 
504 !TEST: version from MOLECULAR VIBRATIONS EB Wilson
505 
506  n23 = sqrt(rpt32(1)*rpt32(1)+rpt32(2)*rpt32(2)+rpt32(3)*rpt32(3))
507  rpt(:) = cp1232(:)*n23/n1/n1
508 !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3))  > tol10) then
509 !write(std_out,*) 'Compare bb1 : '
510 !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
511 !end if
512  bb(:) = rpt(:)
513 
514 end subroutine ddihedral_d1

ABINIT/ddihedral_d2 [ Functions ]

[ Top ] [ Functions ]

NAME

 ddihedral_d2

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

533 subroutine ddihedral_d2(r1,r2,r3,r4,bb)
534 
535  use defs_basis
536  use m_profiling_abi
537 
538  use m_geometry,  only : acrossb
539 
540 !This section has been created automatically by the script Abilint (TD).
541 !Do not modify the following lines by hand.
542 #undef ABI_FUNC
543 #define ABI_FUNC 'ddihedral_d2'
544 !End of the abilint section
545 
546  implicit none
547 
548 !Arguments ------------------------------------
549 !arrays
550  real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3)
551  real(dp),intent(out) :: bb(3)
552 
553 !Local variables
554 !scalars
555  real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,sp1232,sp3432,tmp
556 !arrays
557  real(dp) :: cp1232(3),cp1232_12(3),cp1232_34(3),cp32_1232(3),cp32_3432(3)
558  real(dp) :: cp3432(3),cp3432_12(3),cp3432_34(3),cpcp(3),rpt(3),rpt12(3)
559  real(dp) :: rpt32(3),rpt34(3)
560 
561 ! *************************************************************************
562  rpt12(:) = r1(:)-r2(:)
563  rpt32(:) = r3(:)-r2(:)
564  rpt34(:) = r3(:)-r4(:)
565 
566  call acrossb(rpt12,rpt32,cp1232)
567  call acrossb(rpt34,rpt32,cp3432)
568 
569 !DEBUG
570 !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232
571 !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432
572 !ENDDEBUG
573 
574  n1 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
575  n2 = sqrt(cp3432(1)**2+cp3432(2)**2+cp3432(3)**2)
576 
577  cos_dihedral = (cp1232(1)*cp3432(1)+cp1232(2)*cp3432(2)+cp1232(3)*cp3432(3))/n1/n2
578  if (cos_dihedral > one - epsilon(one)*two) then
579    cos_dihedral = one
580  else if(cos_dihedral < -one + epsilon(one)*two) then
581    cos_dihedral = -one
582  end if
583 !we use complementary of standard angle, so
584 !cos_dihedral = -cos_dihedral
585 
586  call acrossb(cp1232,cp3432,cpcp)
587  cpcp(:) = cpcp(:)/n1/n2
588 !we use complementary of standard angle, but sin is invariant
589  sin_dihedral = -(cpcp(1)*rpt32(1)+cpcp(2)*rpt32(2)+cpcp(3)*rpt32(3))&
590 & /sqrt(rpt32(1)**2+rpt32(2)**2+rpt32(3)**2)
591  dih_sign = one
592  if (sin_dihedral <  -tol12) then
593    dih_sign = -one
594  end if
595 
596 !DEBUG
597 !write(std_out,'(a,3E16.6)') 'ddihedral_d2 : cos abs(sin) dih_sign= ',&
598 !&    cos_dihedral,sin_dihedral,dih_sign
599 !ENDDEBUG
600 
601 !ddihedral_d2 = dih_sign* acos(cos_dihedral)
602  call acrossb(rpt32,cp3432,cp32_3432)
603  call acrossb(cp3432,rpt12,cp3432_12)
604  call acrossb(cp1232,rpt34,cp1232_34)
605 
606  call acrossb(rpt32,cp1232,cp32_1232)
607  call acrossb(cp1232,rpt12,cp1232_12)
608  call acrossb(cp3432,rpt34,cp3432_34)
609 
610  rpt(:) = -(cp32_3432(:) + cp3432_12(:) + cp1232_34(:))/n1/n2 &
611 & +cos_dihedral*(cp32_1232(:)/n1/n1 + cp1232_12(:)/n1/n1 + cp3432_34(:)/n2/n2)
612  bb(:) = zero
613  tmp = sqrt(one-cos_dihedral**2)
614  if (tmp > tol12) then
615 !  we use complementary of standard angle, so cosine in acos has - sign,
616 !  and it appears for derivative
617    bb(:) = -dih_sign * rpt(:) * (-one) / tmp
618  else
619    bb(:) = dih_sign * cos_dihedral * &
620 &   ( cp32_1232(:)/n1/n1/sqrt(cp32_1232(1)**2+cp32_1232(2)**2+cp32_1232(3)**2) &
621 &   +cp1232_12(:)/n1/n1/sqrt(cp1232_12(1)**2+cp1232_12(2)**2+cp1232_12(3)**2) &
622 &   +cp3432_34(:)/n2/n2/sqrt(cp3432_34(1)**2+cp3432_34(2)**2+cp3432_34(3)**2) )
623  end if
624 
625 !TEST: version from MOLECULAR VIBRATIONS EB Wilson p. 61
626  n23 = sqrt(rpt32(1)*rpt32(1)+rpt32(2)*rpt32(2)+rpt32(3)*rpt32(3))
627  sp1232 = rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3)
628  sp3432 = rpt34(1)*rpt32(1)+rpt34(2)*rpt32(2)+rpt34(3)*rpt32(3)
629 
630  rpt(:) = -cp1232(:)*(n23-sp1232/n23)/n1/n1 - cp3432(:)*sp3432/n23/n2/n2
631 !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3))  > tol10) then
632 !write(std_out,*) 'Compare bb2 : '
633 !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
634 !write(std_out,*) -cp1232(:)*(n23-sp1232/n23)/n1/n1, -cp3432(:)*sp3432/n23/n2/n2
635 !end if
636  bb(:) = rpt(:)
637 
638 
639 end subroutine ddihedral_d2