TABLE OF CONTENTS

ABINIT/glue_lotf [ Modules ]

[ Top ] [ Modules ]

NAME

 glue_lotf

FUNCTION

  Contains the GLUE procedure and parameters for Lotf

COPYRIGHT

  Copyright (C) 2009-2024 ABINIT group (MMancini)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  Parameters for the glue potential
  ref: Ercolessi et al, Phil Mag A 58(1), 213 (1988)

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module glue_lotf
27 
28  use defs_basis
29  use m_errors
30 
31  implicit none
32  private
33 
34  !--Parameters for the glue potential (Ercolessi et al, Phil Mag A 58(1), 213 (1988))
35 
36  !--Function phi_r
37  real(dp),public :: dphi,rcphi
38  real(dp) :: a0I,a1I,a2I,a3I,a4I
39  real(dp) :: a0II,a1II,a2II,a3II,a4II,a5II,a6II
40 
41  !--Function rho_r
42  real(dp),public :: rmrho
43  real(dp) :: rbrho
44  real(dp) :: b0I,b1I,b2I,b3I
45  real(dp) :: b0II,b1II,b2II,b3II,b0III,b1III,b2III,b3III
46 
47  !--Function U_n
48  real(dp) :: n0U,nsU
49  real(dp) :: c0I,c1I,c2I,c3I,c4I
50  real(dp) :: c0II,c1II,c2II,c3II,c4II,c0III,c1III,c2III,c3III
51 
52  !--Auxiliary variables to change glue parameters
53  real(dp),parameter :: ddphi = 0.028779246_dp
54  real(dp),parameter :: scalefactor_phi = one
55 
56 
57  public ::             &
58    glue_init,          &
59    glue_pair_devs,     &
60    glue_pair,          &
61    calc_coord,         &
62    calc_rhop,          &
63    calc_coord_new_d,   &
64    rho_devs,           &
65    rhop_value,         &
66    eval_U_n,           &
67    eval_Upp_n
68 
69 contains !===========================================================

eval_lotf/eval_Upp_n [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 eval_Upp_n

FUNCTION

INPUTS

SOURCE

582  subroutine eval_Upp_n(coordatom_i,up_dum,upp_dum)
583 
584   implicit none
585 
586   !Arguments ------------------------
587   real(dp),intent(in) :: coordatom_i
588   real(dp),intent(out) :: up_dum,upp_dum
589   !Local ---------------------------
590   real(dp) :: cns,cns2,cns3,cn0,cn02,cn03
591   character(len=500) :: msg
592 
593 ! *************************************************************************
594 
595   !--Evaluate U and Up for atom i
596   up_dum = zero
597   upp_dum = zero
598 
599   if (coordatom_i < zero) then
600 
601     write(msg,'(3a,i8,2a)')&
602       &    'ERROR: negative coordination!!! ',ch10,&
603       &    'coord = ', coordatom_i,ch10,&
604       &    'The coordination cannot be negative!'
605     ABI_ERROR(msg)
606 
607   elseif (coordatom_i < nsU) then
608 
609     cns = coordatom_i - nsU
610     cns2 = cns * cns
611     cns3 = cns2 * cns
612 
613     up_dum = 4.d0*c4I*cns3 + 3.d0*c3I*cns2 + 2.d0*c2I*cns + c1I
614     upp_dum = 12.d0*c4I*cns2 + 6.d0*c3I*cns + 2.d0*c2I
615 
616   elseif (coordatom_i < n0U) then
617 
618     cn0 = coordatom_i - n0U
619     cn02 = cn0 * cn0
620     cn03 = cn02 * cn0
621 
622     up_dum = 4.d0*c4II*cn03 + 3.d0*c3II*cn02 + 2.d0*c2II*cn0 + c1II
623     upp_dum = 12.d0*c4II*cn02 + 6.d0*c3II*cn0 + 2.d0*c2II
624 
625   elseif (coordatom_i >= n0U) then
626 
627     cn0 = coordatom_i - n0U
628     cn02 = cn0 * cn0
629 
630     up_dum = 3.d0*c3III*cn02 + 2.d0*c2III*cn0 + c1III
631     upp_dum = 6.d0*c3III*cn0 + 2.d0*c2III
632 
633   endif
634  end subroutine eval_Upp_n
635 
636 end module glue_lotf

glue_lotf/calc_coord [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 calc_coord

FUNCTION

INPUTS

SOURCE

304  subroutine calc_coord(r_au,coordatom_dum)
305 
306   implicit none
307 
308   !Arguments ------------------------
309   real(dp) :: r_au,coordatom_dum
310   !Local ---------------------------
311   real(dp) :: rst,r,rho_neigh
312   ! *********************************************************************
313 !    --Adds the density of a single given neighbour
314 
315      rst = sqrt(r_au)
316      rho_neigh = zero
317 
318      if (rst < dphi) then
319        r = rst - dphi
320        rho_neigh = b3I*r*r*r + b2I*r*r + b1I*r + b0I
321      elseif (rst < rbrho) then
322        r = rst - dphi
323        rho_neigh = b3II*r*r*r + b2II*r*r + b1II*r + b0II
324      elseif (rst < rmrho) then
325        r = rst - rmrho
326        rho_neigh = b3III*r*r*r + b2III*r*r + b1III*r + b0III
327      end if
328 
329      coordatom_dum = coordatom_dum + rho_neigh
330 
331    end subroutine calc_coord

glue_lotf/calc_coord_new_d [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 calc_coord_new_d

FUNCTION

INPUTS

SOURCE

386  subroutine calc_coord_new_d(r_au,alpha_d,coordatom_dum)
387 
388   implicit none
389 
390   !Arguments ------------------------
391   real(dp),intent(in) :: r_au,alpha_d
392   real(dp),intent(out) ::coordatom_dum
393   !Local ---------------------------
394   real(dp) :: rst,r,rho_neigh
395   ! *********************************************************************
396 
397 !    --Adds the density of a single given neighbour
398      rst = sqrt(r_au)
399 
400      rho_neigh = zero
401 
402      if (rst < alpha_d) then
403        r = rst - alpha_d
404        rho_neigh = b3I*r*r*r + b2I*r*r + b1I*r + b0I
405      elseif (rst < rbrho) then
406        r = rst - alpha_d
407        rho_neigh = b3II*r*r*r + b2II*r*r + b1II*r + b0II
408 
409 !      --Restriction on rmrho to keep rho and rhop continuous
410      elseif (rst < (rmrho+alpha_d-dphi)) then
411        r = rst - (rmrho+alpha_d-dphi)
412        rho_neigh = b3III*r*r*r + b2III*r*r + b1III*r + b0III
413      end if
414 
415      coordatom_dum = coordatom_dum + rho_neigh
416    end subroutine calc_coord_new_d

glue_lotf/calc_rhop [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 calc_rhop

FUNCTION

INPUTS

SOURCE

346      subroutine calc_rhop(r_st,rhop_dum)
347 
348      implicit none
349 
350 !    Arguments ------------------------
351      real(dp),intent(in) :: r_st
352      real(dp),intent(out) :: rhop_dum
353 !    Local ---------------------------
354      real(dp) :: r,r2
355 !    *********************************************************************
356 
357      rhop_dum = zero
358 
359      if (r_st < dphi) then
360        r = r_st - dphi
361        r2 = r*r
362        rhop_dum = (3.d0*b3I*r2 + 2.d0*b2I*r + b1I)
363      elseif (r_st < rbrho) then
364        r = r_st - dphi
365        r2 = r*r
366        rhop_dum = (3.d0*b3II*r2 + 2.d0*b2II*r + b1II)
367      elseif (r_st < rmrho) then
368        r = r_st - rmrho
369        r2 = r*r
370        rhop_dum = (3.d0*b3III*r2 + 2.d0*b2III*r + b1III)
371      end if
372    end subroutine calc_rhop

glue_lotf/eval_U_n [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 eval_U_n

FUNCTION

INPUTS

SOURCE

516  subroutine eval_U_n(coordatom_i,epot_dum2,up_dum)
517 
518   implicit none
519 
520   !Arguments ------------------------
521   real(dp),intent(in) :: coordatom_i
522   real(dp),intent(out) :: epot_dum2,up_dum
523   !Local ---------------------------
524   real(dp) :: cns,cns2,cns3,cns4,cn0,cn02,cn03,cn04
525   character(len=500) :: msg
526 
527   !--Evaluate U and Up for atom i
528   epot_dum2 = zero
529   up_dum = zero
530 
531   if (coordatom_i < zero) then
532 
533     write(msg,'(3a,i8,2a)')&
534       &    'ERROR: negative coordination!!! ',ch10,&
535       &    'coord = ', coordatom_i,ch10,&
536       &    'The coordination cannot be negative!'
537     ABI_ERROR(msg)
538 
539   elseif (coordatom_i < nsU) then
540 
541     cns = coordatom_i - nsU
542     cns2 = cns * cns
543     cns3 = cns2 * cns
544     cns4 = cns3 * cns
545 
546     epot_dum2 = c4I*cns4 + c3I*cns3 + c2I*cns2 + c1I*cns + c0I
547     up_dum = 4.d0*c4I*cns3 + 3.d0*c3I*cns2 + 2.d0*c2I*cns + c1I
548 
549   elseif (coordatom_i < n0U) then
550 
551     cn0 = coordatom_i - n0U
552     cn02 = cn0 * cn0
553     cn03 = cn02 * cn0
554     cn04 = cn03 * cn0
555 
556     epot_dum2 = c4II*cn04 + c3II*cn03 + c2II*cn02 + c1II*cn0 + c0II
557     up_dum = 4.d0*c4II*cn03 + 3.d0*c3II*cn02 + 2.d0*c2II*cn0 + c1II
558 
559   elseif (coordatom_i >= n0U) then
560 
561     cn0 = coordatom_i - n0U
562     cn02 = cn0 * cn0
563     cn03 = cn02 * cn0
564 
565     epot_dum2 = c3III*cn03 + c2III*cn02 + c1III*cn0 + c0III
566     up_dum = 3.d0*c3III*cn02 + 2.d0*c2III*cn0 + c1III
567 
568   endif
569  end subroutine eval_U_n

glue_lotf/glue_init [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 glue_init

FUNCTION

INPUTS

SOURCE

 81  subroutine glue_init()
 82 
 83   !--Function phi_r
 84 ! *************************************************************************
 85   dphi =  (0.2878207442141723d+01 + ddphi) / 0.529177d0
 86   rcphi = (0.3700000000000000d+01 + ddphi) / 0.529177d0
 87 
 88   a0I = -0.8000000000000000d-01 / (13.6057d0 * 2.d0) *&
 89     scalefactor_phi
 90   a1I =  0.0000000000000000d+00 / (13.6057d0 * 2.d0) * (0.529177d0) *&
 91     scalefactor_phi
 92   a2I =  0.7619231375231362d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**2 *&
 93     scalefactor_phi
 94   a3I = -0.8333333333333333d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**3 *&
 95     scalefactor_phi
 96   a4I = -0.1211483464993159d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**4 *&
 97     scalefactor_phi
 98 
 99   !        a0II = -0.8000000000000000d-01 / (13.6057d0 * 2.d0)
100   a0II = a0I
101   !        a1II =  0.0000000000000000d+00 / (13.6057d0 * 2.d0) * (0.529177d0)
102   a1II = a1I
103   !        a2II =  0.7619231375231362d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**2
104   a2II = a2I
105   a3II = -0.8333333333333333d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**3 *&
106     scalefactor_phi
107   a4II = -0.1096009851140349d+01 / (13.6057d0 * 2.d0) * (0.529177d0)**4 *&
108     scalefactor_phi
109   a5II =  0.2158417178555998d+01 / (13.6057d0 * 2.d0) * (0.529177d0)**5 *&
110     scalefactor_phi
111   a6II = -0.9128915709636862d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**6 *&
112     scalefactor_phi
113 
114   !--Function rho_r
115   rbrho = 0.3500000000000000d+01 / 0.529177d0
116   rmrho = (0.3900000000000000d+01 + ddphi) / 0.529177d0
117 
118   b0I =  0.1000000000000000d+01
119   b1I = -0.6800000000000000d+00 * (0.529177d0)
120   b2I =  0.7500000000000000d+00 * (0.529177d0)**2
121   b3I = -0.1333333333333333d+01 * (0.529177d0)**3
122 
123   !        b0II =  0.1000000000000000d+01
124   b0II = b0I
125   !        b1II = -0.6800000000000000d+00 * (0.529177d0)
126   b1II = b1I
127   !        b2II =  0.7500000000000000d+00 * (0.529177d0)**2
128   b2II = b2I
129   b3II = -0.1527241171296038d+01 * (0.529177d0)**3
130 
131   b0III =  0.0000000000000000d+00
132   b1III =  0.0000000000000000d+00 * (0.529177d0)
133   b2III =  0.5578188675490974d+01 * (0.529177d0)**2
134   b3III =  0.6132971688727435d+01 * (0.529177d0)**3
135 
136   !--Function U_n
137   n0U =  0.1200000000000000d+02
138   nsU =  0.9358157767784574d+01
139 
140   c0I = -0.2793388616771698d+01 / (13.6057d0 * 2.d0) * scalefactor_phi
141   c1I = -0.3419999999999999d+00 / (13.6057d0 * 2.d0) * scalefactor_phi
142   c2I =  0.3902327808424106d-01 / (13.6057d0 * 2.d0) * scalefactor_phi
143   c3I =  0.7558829951858879d-02 / (13.6057d0 * 2.d0) * scalefactor_phi
144   c4I =  0.3090472511796849d-03 / (13.6057d0 * 2.d0) * scalefactor_phi
145 
146   c0II = -0.3300000000000000d+01 / (13.6057d0 * 2.d0) * scalefactor_phi
147   c1II =  0.0000000000000000d+00 / (13.6057d0 * 2.d0) * scalefactor_phi
148   c2II =  0.8618226772941980d-01 / (13.6057d0 * 2.d0) * scalefactor_phi
149   c3II =  0.4341701445034724d-02 / (13.6057d0 * 2.d0) * scalefactor_phi
150   c4II = -0.3044398779375916d-03 / (13.6057d0 * 2.d0) * scalefactor_phi
151 
152   !        c0III = -0.3300000000000000d+01 / (13.6057d0 * 2.d0)
153   c0III = c0II
154   !        c1III =  0.0000000000000000d+00 / (13.6057d0 * 2.d0)
155   c1III = c1II
156   !        c2III =  0.8618226772941980d-01 / (13.6057d0 * 2.d0)
157   c2III = c2II
158   c3III =  0.4325981467602070d-02 / (13.6057d0 * 2.d0) * scalefactor_phi
159 
160  end subroutine glue_init

glue_lotf/glue_pair [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 glue_pair

FUNCTION

INPUTS

SOURCE

248  subroutine glue_pair(RD,r_au,epot_2,fdum)
249 
250   implicit none
251 
252   !Arguments ------------------------
253   real(dp) ::  RD(3), fdum(3), r_au, epot_2
254   !Local ---------------------------
255   real(dp) ::  rst, r, r2, r3, r4, r5, r6
256   integer ::  k
257   ! *********************************************************************
258 !    --The calculation for a given pair
259      rst = sqrt(r_au)
260      r = rst - dphi
261      r2 = r*r
262      r3 = r2*r
263      r4 = r3*r
264 
265 !    --Energy
266      epot_2 = zero
267      if (rst < dphi) then
268        epot_2 = a4I*r4 + a3I*r3 + a2I*r2 + a1I*r + a0I
269      elseif (rst < rcphi) then
270        r5 = r4*r
271        r6 = r5*r
272        epot_2 = a6II*r6 + a5II*r5 + a4II*r4 + a3II*r3 + a2II*r2 +&
273        a1II*r + a0II
274      end if
275 
276 !    --Forces (without the minus sign, because we use rji)
277      fdum(:) = zero
278 
279      do k =1,3
280        if (rst < dphi) then
281          fdum(k) =  ( 4.d0*a4I*r3 + 3.d0*a3I*r2 + 2.d0*a2I*r + a1I ) * &
282          RD(k)/rst
283        elseif (rst < rcphi) then
284          fdum(k) =  ( 6.d0*a6II*r5 + 5.d0*a5II*r4 + 4.d0*a4II*r3 + &
285          3.d0*a3II*r2 + 2.d0*a2II*r + a1II ) * &
286          RD(k)/rst
287        end if
288      end do
289 
290    end subroutine glue_pair

glue_lotf/glue_pair_devs [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 glue_pair_devs

FUNCTION

INPUTS

SOURCE

173  subroutine glue_pair_devs(alpha_dum,RD,r_au,epot_2,fdum,dfdum)
174 
175   implicit none
176 
177   !Arguments ------------------------
178   real(dp),intent(in) ::  alpha_dum(3),RD(3)
179   real(dp),intent(in) ::  r_au
180   real(dp),intent(out) ::  epot_2
181   real(dp),intent(out) ::  fdum(3),dfdum(3,2)
182   !Local ---------------------------
183   real(dp) :: rst, r, r2, r3, r4, r5, r6
184   real(dp) :: drcphi, dphi2, rcphi2, scale_factor
185   ! *********************************************************************
186 
187 !    --The calculation for a given pair
188      drcphi = rcphi - dphi
189      dphi2 = alpha_dum(1)
190 
191 !    --Rewritten
192      rcphi2 = dphi2 + drcphi
193 
194      scale_factor = alpha_dum(2)
195 
196      rst = sqrt(r_au)
197 
198      r = rst - dphi2
199      r2 = r*r
200      r3 = r2*r
201      r4 = r3*r
202 
203 !    --Energy
204      epot_2 = zero
205 
206      if (rst < dphi2) then
207        epot_2 = scale_factor * (a4I*r4 + a3I*r3 + a2I*r2 + a1I*r + a0I)
208      elseif (rst < rcphi2) then
209        r5 = r4*r
210        r6 = r5*r
211        epot_2 = scale_factor * (a6II*r6 + a5II*r5 + a4II*r4 + a3II*r3 + a2II*r2 +&
212        a1II*r + a0II)
213      end if
214 
215 !    --Forces (without the minus sign, because we use rji) and derivatives w.r.t. parameter d
216      fdum(:) = zero
217      dfdum(:,:) = zero
218 
219      if (rst < dphi2) then
220        fdum(:) = scale_factor * (4.d0*a4I*r3 + 3.d0*a3I*r2 + 2.d0*a2I*r + a1I) * RD(:)/rst
221        dfdum(:,1) = -scale_factor * (12.d0*a4I*r2 + 6.d0*a3I*r + 2.d0*a2I) * RD(:)/rst
222        dfdum(:,2) = fdum(:) / scale_factor
223 
224      elseif (rst < rcphi2) then
225        fdum(:) = (scale_factor) * &
226 &       (6.d0*a6II*r5 + 5.d0*a5II*r4 + 4.d0*a4II*r3 + 3.d0*a3II*r2 + 2.d0*a2II*r + a1II) * &
227 &       RD(:)/rst
228 
229        dfdum(:,1) = -scale_factor * &
230 &       (30.d0*a6II*r4 + 20.d0*a5II*r3 + 12.d0*a4II*r2 + 6.d0*a3II*r + 2.d0*a2II) * RD(:)/rst
231 
232        dfdum(:,2) = fdum(:) / scale_factor
233 
234      end if
235    end subroutine glue_pair_devs

glue_lotf/rho_devs [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 rho_devs

FUNCTION

INPUTS

SOURCE

429      subroutine rho_devs(r_au,alpha_d,rho_neigh_d,rho_neigh_p,rho_neigh_pd)
430 
431      real(dp),intent(in) :: r_au,alpha_d
432      real(dp),intent(out) :: rho_neigh_d,rho_neigh_p,rho_neigh_pd
433 !    Local ---------------------------
434      real(dp) :: rst,r
435 !    *********************************************************************
436 
437 !    --Adds the density of a single given neighbour
438      rst = sqrt(r_au)
439 
440      rho_neigh_p = zero
441      rho_neigh_d = zero
442      rho_neigh_pd = zero
443 
444      if (rst < alpha_d) then
445        r = rst - alpha_d
446        rho_neigh_d = - ( 3.d0*b3I*r*r + 2.d0*b2I*r + b1I )
447        rho_neigh_p = - rho_neigh_d
448        rho_neigh_pd = - ( 6.d0*b3I*r + 2.d0*b2I )
449 
450      elseif (rst < rbrho) then
451        r = rst - alpha_d
452        rho_neigh_d = - ( 3.d0*b3II*r*r + 2.d0*b2II*r + b1II )
453        rho_neigh_p = - rho_neigh_d
454        rho_neigh_pd = - ( 6.d0*b3II*r + 2.d0*b2II )
455 
456 !      --Restriction on rmrho to keep rho and rhop continuous
457      elseif (rst < (rmrho+alpha_d-dphi)) then
458        r = rst - (rmrho+alpha_d-dphi)
459        rho_neigh_d = zero
460        rho_neigh_p = 3.d0*b3III*r*r + 2.d0*b2III*r + b1III
461        rho_neigh_pd = zero
462 
463      end if
464    end subroutine rho_devs

glue_lotf/rhop_value [ Functions ]

[ Top ] [ glue_lotf ] [ Functions ]

NAME

 rhop_value

FUNCTION

INPUTS

SOURCE

478  subroutine rhop_value(rst,alpha_d,rhop_dum)
479 
480   implicit none
481 
482   !Arguments ------------------------
483   real(dp),intent(in) :: rst,alpha_d
484   real(dp),intent(out) :: rhop_dum
485   !Local ---------------------------
486   real(dp) :: r
487   ! *********************************************************************
488 
489 !    --Adds the density of a single given neighbour
490      rhop_dum = zero
491 
492      if (rst < alpha_d) then
493        r = rst - alpha_d
494        rhop_dum = 3.d0*b3I*r*r + 2.d0*b2I*r + b1I
495      elseif (rst < rbrho) then
496        r = rst - alpha_d
497        rhop_dum = 3.d0*b3II*r*r + 2.d0*b2II*r + b1II
498 !      --Restriction on rmrho to keep rho and rhop continuous
499      elseif (rst < (rmrho+alpha_d-dphi)) then
500        r = rst - (rmrho+alpha_d-dphi)
501        rhop_dum = 3.d0*b3III*r*r + 2.d0*b2III*r + b1III
502      end if
503    end subroutine rhop_value