RMPCDMD
md.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2015-2017 Pierre de Buyl and contributors
3 ! License: BSD 3-clause (see file LICENSE)
4 
12 
13 module md
14  use particle_system
15  use cell_system
16  use interaction
17  use common
18  use quaternion
19  implicit none
20 
21  private
22 
23  public :: md_pos
24  public :: apply_pbc
25  public :: md_vel
26  public :: rattle_dimer_pos
27  public :: rattle_dimer_vel
29  public :: lj93_zwall
30  public :: elastic_network
31  public :: rigid_body_t
32  public :: cell_md_pos, cell_md_vel
33  public :: cell_md_pos_ywall
34  public :: cell_md_pos_zwall
35 
40  integer :: i_start
41  integer :: i_stop
42  double precision :: mass
43  double precision :: pos(3)
44  double precision, allocatable :: pos_body(:,:)
45  double precision :: vel(3)
46  double precision :: force(3)
47  double precision :: q(4)
48  double precision :: l(3)
49  double precision :: torque(3)
50  double precision :: i_body(3)
51  double precision :: omega_body(3)
52  contains
53  procedure :: init => rigid_body_init
54  procedure :: compute_force_torque => rigid_body_compute_force_torque
55  procedure :: vv1 => rigid_body_vv1
56  procedure :: vv2 => rigid_body_vv2
57  end type rigid_body_t
58 
59 contains
60 
61  subroutine md_pos(particles, dt)
62  type(particle_system_t), intent(inout) :: particles
63  double precision, intent(in) :: dt
64 
65  double precision :: dt_sq
66  integer :: k
67 
68  dt_sq = dt**2/2
69 
70  call particles%time_md_pos%tic()
71  !$omp parallel do
72  do k = 1, particles% Nmax
73  particles% pos(:,k) = particles% pos(:,k) + dt * particles% vel(:,k) + dt_sq * particles% force(:,k)
74  end do
75  call particles%time_md_pos%tac()
76 
77  end subroutine md_pos
78 
79  subroutine cell_md_pos(cells, particles, dt, md_flag)
80  type(cell_system_t), intent(inout) :: cells
81  type(particle_system_t), intent(inout) :: particles
82  double precision, intent(in) :: dt
83  logical, intent(in) :: md_flag
84 
85  double precision :: dt_sq
86  integer :: i, j
87 
88  dt_sq = dt**2/2
89 
90  call particles%time_md_pos%tic()
91  !$omp parallel do private(i, j)
92  do i = 1, cells%N
93  if (cells%is_md(i) .eqv. md_flag) then
94  do j = cells%cell_start(i), cells%cell_start(i) + cells%cell_count(i) - 1
95  particles% pos(:,j) = particles% pos(:,j) + dt * particles% vel(:,j) + dt_sq * particles% force(:,j)
96  end do
97  end if
98  end do
99 
100  call particles%time_md_pos%tac()
101 
102  end subroutine cell_md_pos
103 
104  subroutine cell_md_pos_ywall(cells, particles, dt, md_flag)
105  type(cell_system_t), intent(inout) :: cells
106  type(particle_system_t), intent(inout) :: particles
107  double precision, intent(in) :: dt
108  logical, intent(in) :: md_flag
109 
110  double precision :: dt_sq
111  double precision :: Ly, x(3), v(3), overtime
112  integer :: i, j, bc
113  logical :: collide
114 
115  dt_sq = dt**2/2
116  ly = cells%edges(2)
117  bc = cells%bc(2)
118 
119  call particles%time_md_pos%tic()
120  !$omp parallel do private(i, j, x, v, collide, overtime)
121  do i = 1, cells%N
122  if (cells%is_md(i) .eqv. md_flag) then
123  do j = cells%cell_start(i), cells%cell_start(i) + cells%cell_count(i) - 1
124  v = particles%vel(:,j)
125  x = particles%pos(:,j) + dt * v + dt_sq * particles%force(:,j)
126  collide = .false.
127  if (x(2)<0.d0) then
128  collide = .true.
129  overtime = x(2)/v(2)
130  else if (x(2)>ly) then
131  collide = .true.
132  overtime = (x(2)-ly)/v(2)
133  end if
134  if (collide) then
135  if (bc==bounce_back_bc) then
136  v = -v
137  x = x+2*overtime*v
138  else if (bc==specular_bc) then
139  v(2) = -v(2)
140  x(2) = x(2)+2*overtime*v(2)
141  end if
142  particles%vel(:,j) = v
143  end if
144  particles% pos(:,j) = x
145  end do
146  end if
147  end do
148 
149  call particles%time_md_pos%tac()
150 
151  end subroutine cell_md_pos_ywall
152 
153  subroutine cell_md_pos_zwall(cells, particles, dt, md_flag)
154  type(cell_system_t), intent(inout) :: cells
155  type(particle_system_t), intent(inout) :: particles
156  double precision, intent(in) :: dt
157  logical, intent(in) :: md_flag
158 
159  double precision :: dt_sq
160  double precision :: Lz, x(3), v(3), overtime
161  integer :: i, j, bc
162  logical :: collide
163 
164  dt_sq = dt**2/2
165  lz = cells%edges(3)
166  bc = cells%bc(3)
167 
168  call particles%time_md_pos%tic()
169  !$omp parallel do private(i, j, x, v, collide, overtime)
170  do i = 1, cells%N
171  if (cells%is_md(i) .eqv. md_flag) then
172  do j = cells%cell_start(i), cells%cell_start(i) + cells%cell_count(i) - 1
173  v = particles%vel(:,j)
174  x = particles%pos(:,j) + dt * v + dt_sq * particles%force(:,j)
175  collide = .false.
176  if (x(3)<0.d0) then
177  collide = .true.
178  overtime = x(3)/v(3)
179  else if (x(3)>lz) then
180  collide = .true.
181  overtime = (x(3)-lz)/v(3)
182  end if
183  if (collide) then
184  if (bc==bounce_back_bc) then
185  v = -v
186  x = x+2*overtime*v
187  else if (bc==specular_bc) then
188  v(3) = -v(3)
189  x(3) = x(3)+2*overtime*v(3)
190  end if
191  particles%vel(:,j) = v
192  end if
193  particles% pos(:,j) = x
194  end do
195  end if
196  end do
197 
198  call particles%time_md_pos%tac()
199 
200  end subroutine cell_md_pos_zwall
201 
202  subroutine apply_pbc(particles, edges)
203  type(particle_system_t), intent(inout) :: particles
204  double precision, intent(in) :: edges(3)
205 
206  integer :: k, jump(3)
207 
208  call particles%time_apply_pbc%tic()
209  !$omp parallel do private(jump)
210  do k = 1, particles% Nmax
211  jump = floor(particles% pos(:,k) / edges)
212  particles% image(:,k) = particles% image(:,k) + jump
213  particles% pos(:,k) = particles% pos(:,k) - jump*edges
214  end do
215  call particles%time_apply_pbc%tac()
216 
217  end subroutine apply_pbc
218 
219  subroutine md_vel(particles, dt)
220  type(particle_system_t), intent(inout) :: particles
221  double precision, intent(in) :: dt
222 
223  integer :: k
224 
225  call particles%time_md_vel%tic()
226  !$omp parallel do
227  do k = 1, particles% Nmax
228  if (.not. btest(particles%flags(k), wall_bit)) then
229  particles% vel(:,k) = particles% vel(:,k) + &
230  dt * ( particles% force(:,k) + particles% force_old(:,k) ) / 2
231  else
232  particles%flags(k) = ibclr(particles%flags(k), wall_bit)
233  end if
234  end do
235  call particles%time_md_vel%tac()
236 
237  end subroutine md_vel
238 
239  subroutine cell_md_vel(cells, particles, dt, md_flag)
240  type(cell_system_t), intent(in) :: cells
241  type(particle_system_t), intent(inout) :: particles
242  double precision, intent(in) :: dt
243  logical, intent(in) :: md_flag
244 
245  integer :: i, j
246 
247  call particles%time_md_vel%tic()
248  !$omp parallel do private(i, j)
249  do i = 1, cells%N
250  if (cells%is_md(i) .eqv. md_flag) then
251  do j = cells%cell_start(i), cells%cell_start(i) + cells%cell_count(i) - 1
252  if (.not. btest(particles%flags(j), wall_bit)) then
253  particles% vel(:,j) = particles% vel(:,j) + &
254  dt * ( particles% force(:,j) + particles% force_old(:,j) ) / 2
255  else
256  particles%flags(j) = ibclr(particles%flags(j), wall_bit)
257  end if
258  end do
259  end if
260  end do
261  call particles%time_md_vel%tac()
262 
263  end subroutine cell_md_vel
264 
265  subroutine rattle_dimer_pos(p, d, dt,edges)
266  type(particle_system_t), intent(inout) :: p
267  double precision, intent(in) :: d
268  double precision, intent(in) :: dt
269  double precision, intent(in) :: edges(3)
270 
271  double precision :: g
272  double precision :: s(3) ! direction vector
273  double precision :: r(3) ! old direction vector
274  double precision :: rsq, ssq, rs
275  double precision :: mass1, mass2, inv_mass
276 
277  r = rel_pos(p% pos_rattle(:,1),p% pos_rattle(:,2), edges)
278  s = rel_pos(p% pos(:,1),p% pos(:,2), edges)
279  mass1 = p%mass(p%species(1))
280  mass2 = p%mass(p%species(2))
281  inv_mass = 1/mass1 + 1/mass2
282 
283  rsq = dot_product(r,r)
284  ssq = dot_product(s,s)
285  rs = dot_product(r,s)
286 
287  g = rs - sqrt(rs**2 - rsq*(ssq-d**2))
288  g = g / (dt * inv_mass * rsq)
289 
290  p% pos(:,1) = p% pos(:,1) - g*dt*r/mass1
291  p% pos(:,2) = p% pos(:,2) + g*dt*r/mass2
292 
293  p% vel(:,1) = p% vel(:,1) - g*r/mass1
294  p% vel(:,2) = p% vel(:,2) + g*r/mass2
295 
296  end subroutine rattle_dimer_pos
297 
298  subroutine rattle_dimer_vel(p, d, dt,edges)
299  type(particle_system_t), intent(inout) :: p
300  double precision, intent(in) :: d
301  double precision, intent(in) :: dt
302  double precision, intent(in) :: edges(3)
303 
304  double precision :: k !second correction factor
305  double precision :: s(3) !direction vector
306  double precision :: mass1, mass2, inv_mass
307 
308  mass1 = p%mass(p%species(1))
309  mass2 = p%mass(p%species(2))
310  inv_mass = 1/mass1 + 1/mass2
311 
312  s = rel_pos(p% pos(:,1), p% pos(:,2), edges)
313 
314  k = dot_product(p%vel(:,1)-p%vel(:,2), s) / (d**2*inv_mass)
315 
316  p% vel(:,1) = p% vel(:,1) - k*s/mass1
317  p% vel(:,2) = p% vel(:,2) + k*s/mass2
318 
319  end subroutine rattle_dimer_vel
320 
321  subroutine rattle_body_pos(p, links, distances, dt, edges, precision)
322  type(particle_system_t), intent(inout) :: p
323  integer, intent(in) :: links(:,:)
324  double precision, intent(in) :: distances(:)
325  double precision, intent(in) :: dt, edges(3), precision
326 
327  double precision :: g, d, error
328  double precision :: s(3) ! direction vector
329  double precision :: r(3) ! old direction vector
330  double precision :: rsq, ssq, rs
331  double precision :: i_mass1, i_mass2, inv_mass
332  double precision, allocatable :: i_mass(:)
333  integer, allocatable :: im(:,:), im_r(:,:)
334 
335  integer :: rattle_i, rattle_max, i_link, n_link
336  integer :: i1, i2
337 
338  n_link = size(links, dim=2)
339  rattle_max = 1000
340 
341  allocate(i_mass(size(p%mass)))
342  i_mass = 1/p%mass
343 
344  allocate(im(3,p%Nmax))
345  allocate(im_r(3,p%Nmax))
346 
347  s = p%pos(:,1)
348  do i1 = 2, p%Nmax
349  im(:,i1) = floor(((p%pos(:,i1)-s)/edges) + 0.5d0)
350  p%pos(:,i1) = p%pos(:,i1) - im(:,i1)*edges
351  end do
352 
353  s = p%pos_rattle(:,1)
354  do i1 = 2, p%Nmax
355  im_r(:,i1) = floor(((p%pos_rattle(:,i1)-s)/edges) + 0.5d0)
356  p%pos_rattle(:,i1) = p%pos_rattle(:,i1) - im_r(:,i1)*edges
357  end do
358 
359  call p%time_rattle_pos%tic()
360  rattle_loop: do rattle_i = 1, rattle_max
361  error = 0
362  do i_link = 1, n_link
363  i1 = links(1,i_link)
364  i2 = links(2,i_link)
365  d = distances(i_link)
366 
367  r = p%pos_rattle(:,i1)-p%pos_rattle(:,i2)
368  s = p%pos(:,i1)-p%pos(:,i2)
369  i_mass1 = i_mass(p%species(i1))
370  i_mass2 = i_mass(p%species(i2))
371  inv_mass = i_mass1 + i_mass2
372 
373  rsq = dot_product(r,r)
374  ssq = dot_product(s,s)
375  rs = dot_product(r,s)
376 
377  g = rs - sqrt(rs**2 - rsq*(ssq-d**2))
378  g = g / (dt * inv_mass * rsq)
379 
380  p% pos(:,i1) = p% pos(:,i1) - g*dt*r*i_mass1
381  p% pos(:,i2) = p% pos(:,i2) + g*dt*r*i_mass2
382 
383  p% vel(:,i1) = p% vel(:,i1) - g*r*i_mass1
384  p% vel(:,i2) = p% vel(:,i2) + g*r*i_mass2
385  end do
386 
387  do i_link = 1, n_link
388  i1 = links(1,i_link)
389  i2 = links(2,i_link)
390  d = distances(i_link)
391  r = p%pos(:,i1)-p%pos(:,i2)
392  g = abs(sqrt(dot_product(r,r)) - d)
393  if (g > error) error = g
394  end do
395  if (error < precision) exit rattle_loop
396 
397  end do rattle_loop
398  call p%time_rattle_pos%tac()
399  deallocate(i_mass)
400 
401  do i1 = 2, p%Nmax
402  p%pos(:,i1) = p%pos(:,i1) + im(:,i1)*edges
403  end do
404  deallocate(im, im_r)
405 
406  if (rattle_i==rattle_max) write(*,*) 'rattle_max reached in rattle_body_pos'
407 
408  end subroutine rattle_body_pos
409 
410  subroutine rattle_body_vel(p, links, distances, dt, edges, precision)
411  type(particle_system_t), intent(inout) :: p
412  integer, intent(in) :: links(:,:)
413  double precision, intent(in) :: distances(:)
414  double precision, intent(in) :: dt, edges(3), precision
415 
416  double precision :: g, d, error
417  double precision :: s(3), k
418  double precision :: i_mass1, i_mass2, inv_mass
419  double precision, allocatable :: i_mass(:)
420  integer, allocatable :: im(:,:)
421 
422  integer :: rattle_i, rattle_max, i_link, n_link
423  integer :: i1, i2
424 
425  n_link = size(links, dim=2)
426  rattle_max = 1000
427  allocate(i_mass(size(p%mass)))
428  i_mass = 1/p%mass
429 
430  allocate(im(3,p%Nmax))
431 
432  s = p%pos(:,1)
433  do i1 = 2, p%Nmax
434  im(:,i1) = floor(((p%pos(:,i1)-s)/edges) + 0.5d0)
435  p%pos(:,i1) = p%pos(:,i1) - im(:,i1)*edges
436  end do
437 
438  call p%time_rattle_vel%tic()
439  rattle_loop: do rattle_i = 1, rattle_max
440  error = 0
441  do i_link = 1, n_link
442  i1 = links(1,i_link)
443  i2 = links(2,i_link)
444  d = distances(i_link)
445  i_mass1 = i_mass(p%species(i1))
446  i_mass2 = i_mass(p%species(i2))
447  inv_mass = i_mass1 + i_mass2
448 
449  s = p%pos(:,i1)-p%pos(:,i2)
450 
451  k = dot_product(p%vel(:,i1)-p%vel(:,i2), s) / (d**2*inv_mass)
452 
453  p% vel(:,i1) = p% vel(:,i1) - k*s*i_mass1
454  p% vel(:,i2) = p% vel(:,i2) + k*s*i_mass2
455  end do
456 
457  do i_link = 1, n_link
458  i1 = links(1,i_link)
459  i2 = links(2,i_link)
460  d = distances(i_link)
461  i_mass1 = i_mass(p%species(i1))
462  i_mass2 = i_mass(p%species(i2))
463  inv_mass = i_mass1 + i_mass2
464  s = p%pos(:,i1)-p% pos(:,i2)
465  k = abs(dot_product(p%vel(:,i1)-p%vel(:,i2), s) / (d**2*inv_mass))
466  if (k>error) error = k
467  end do
468  if (error < precision) exit rattle_loop
469  end do rattle_loop
470  call p%time_rattle_vel%tac()
471 
472  do i1 = 2, p%Nmax
473  p%pos(:,i1) = p%pos(:,i1) + im(:,i1)*edges
474  end do
475 
476  deallocate(i_mass, im)
477 
478  if (rattle_i==rattle_max) write(*,*) 'rattle_max reached in rattle_body_vel'
479 
480  end subroutine rattle_body_vel
481 
482  function lj93_zwall(particles, edges, lj_params) result(e)
483  type(particle_system_t), intent(inout) :: particles
484  double precision, intent(in) :: edges(3)
485  type(lj_params_t), intent(in) :: lj_params
486  double precision :: e
487 
488  integer :: i, s, dim
489  double precision :: z, z_sq, f, dir
490 
491  e = 0
492  do i = 1, particles%Nmax
493  s = particles%species(i)
494  if (s<=0) cycle
495  do dim = 1, 3
496  if (lj_params%sigma(dim, s)<0) cycle
497  z = particles%pos(dim,i)
498  if (z > edges(dim)/2) then
499  z = edges(dim) - z - lj_params%shift(dim)
500  dir = -1
501  else
502  dir = 1
503  z = z - lj_params%shift(dim)
504  end if
505  z_sq = z**2
506  if (z_sq <= lj_params% cut_sq(dim,s)) then
507  f = lj_force_9_3(z, z_sq, lj_params%epsilon(dim,s), lj_params%sigma(dim,s))
508  particles%force(dim,i) = particles%force(dim,i) + dir*f
509  e = e + lj_energy_9_3(z_sq, lj_params%epsilon(dim,s), lj_params%sigma(dim,s))
510  end if
511  end do
512  end do
513 
514  end function lj93_zwall
515 
516  function elastic_network(p, links, distances, k, edges) result(e)
517  type(particle_system_t), intent(inout) :: p
518  integer, intent(in) :: links(:,:)
519  double precision, intent(in) :: distances(:)
520  double precision, intent(in) :: k
521  double precision, intent(in) :: edges(3)
522  double precision :: e
523 
524  integer :: i, i1, i2, n_links
525  double precision :: x(3), f(3), r, r0
526 
527  n_links = size(distances)
528 
529  call p%time_elastic%tic()
530  e = 0
531  do i = 1, n_links
532  i1 = links(1,i)
533  i2 = links(2,i)
534  r0 = distances(i)
535  x = rel_pos(p%pos(:,i1), p%pos(:,i2), edges)
536  r = norm2(x)
537  e = e + k*(r-r0)**2/2
538  f = - k * (r - r0) * x / r
539  p%force(:,i1) = p%force(:,i1) + f
540  p%force(:,i2) = p%force(:,i2) - f
541  end do
542  call p%time_elastic%tac()
543 
544  end function elastic_network
545 
546  !! Compute body frame positions and inertia tensor
547  subroutine rigid_body_init(this, ps, i_start, i_stop, edges)
548  class(rigid_body_t), intent(out) :: this
549  type(particle_system_t), intent(inout) :: ps
550  integer, intent(in) :: i_start
551  integer, intent(in) :: i_stop
552  double precision, intent(in) :: edges(3)
553 
554  double precision :: mass
555  double precision :: pos(3)
556  integer :: i
557 
558  this%i_start = i_start
559  this%i_stop = i_stop
560 
561  this%mass = 0
562  this%pos = 0
563  this%vel = 0
564  do i = i_start, i_stop
565  mass = ps%mass(ps%species(i))
566  this%mass = this%mass + mass
567  this%pos = this%pos + mass*(ps%pos(:,i)+ps%image(:,i)*edges)
568  this%vel = this%vel + mass*ps%vel(:,i)
569  end do
570  this%pos = this%pos/this%mass
571  this%vel = this%vel/this%mass
572 
573  allocate(this%pos_body(3,i_stop-i_start+1))
574 
575  this%I_body = 0
576  do i = i_start, i_stop
577  mass = ps%mass(ps%species(i))
578  pos = ps%pos(:,i)+ps%image(:,i)*edges - this%pos
579  this%pos_body(:,i) = pos
580  this%I_body = this%I_body + mass*(sum(pos**2)-pos**2)
581  end do
582 
583  this%force = 0
584  this%torque = 0
585  this%q = qnew(s=1.d0)
586  this%L = 0
587 
588  end subroutine rigid_body_init
589 
590  !! Compute torque and force on the rigid body in the lab frame
591  subroutine rigid_body_compute_force_torque(this, ps)
592  class(rigid_body_t), intent(inout) :: this
593  type(particle_system_t), intent(in) :: ps
594 
595  double precision :: pos(3), f(3)
596  integer :: i
597 
598  this%force = 0
599  this%torque = 0
600  do i = this%i_start, this%i_stop
601  f = ps%force(:,i)
602  this%force = this%force + f
603  pos = qrot(this%q, this%pos_body(:,i))
604  this%torque = this%torque + cross(pos, f)
605  end do
606 
607  end subroutine rigid_body_compute_force_torque
608 
613  subroutine rigid_body_vv1(this, ps, dt, treshold, edges)
614  class(rigid_body_t), intent(inout) :: this
615  type(particle_system_t), intent(inout) :: ps
616  double precision, intent(in) :: dt
617  double precision, intent(in) :: treshold
618  double precision, intent(in) :: edges(3)
619 
620  double precision :: L(3), L_body(3), L_body_dot(3), omega_body(3), q(4), q_old(4), q_dot(4), torque_body(3)
621  double precision :: I_inv(3)
622  double precision :: old_pos(3), delta_pos(3)
623 
624  integer :: i, j
625 
626  i_inv = 1/this%I_body
627 
628  ! timestep 0
629  l = this%L
630  q = this%q
631 
632  l_body = qrot(qconj(q), l)
633  torque_body = qrot(qconj(q), this%torque)
634  omega_body = l_body*i_inv
635 
636  l_body_dot = torque_body - cross(omega_body, l_body)
637 
638  ! timestep 1/2
639  l_body = l_body + dt*l_body_dot/2
640  omega_body = l_body*i_inv
641  q_dot = qmul(q, qnew(v=omega_body))/2
642  q = q + q_dot*dt/2
643  q = qnormalize(q)
644  l = this%L + this%torque*dt/2
645 
646  ! Iterative procedure
647 
648  q_old = q
649  iterative_q: do while (.true.)
650  l_body = qrot(qconj(q), l)
651  omega_body = i_inv*l_body
652  q_dot = qmul(q, qnew(v=omega_body))/2
653  q = this%q + q_dot*dt/2
654  q = qnormalize(q)
655  if (qnorm(q-q_old)<treshold) exit iterative_q
656  q_old = q
657  end do iterative_q
658 
659  this%q = qnormalize(this%q + dt*q_dot)
660 
661  this%pos = this%pos + this%vel*dt + this%force*dt**2/(2*this%mass)
662 
663  do i = this%i_start, this%i_stop
664  old_pos = ps%pos(:,i)
665  ps%pos(:,i) = this%pos + qrot(this%q, this%pos_body(:,i))
666  delta_pos = ps%pos(:,i) - old_pos + edges/2
667  ps%pos(:,i) = ps%pos(:,i) - floor(delta_pos / edges)*edges
668  end do
669 
670  this%L = this%L + this%torque*dt/2
671  this%vel = this%vel + this%force*dt/(2*this%mass)
672 
673  end subroutine rigid_body_vv1
674 
679  subroutine rigid_body_vv2(this, ps, dt)
680  class(rigid_body_t), intent(inout) :: this
681  type(particle_system_t), intent(inout) :: ps
682  double precision, intent(in) :: dt
683 
684  double precision :: L_body(3), omega(3)
685  double precision :: pos(3)
686  integer :: i
687 
688  this%L = this%L + this%torque*dt/2
689  l_body = qrot(qconj(this%q), this%L)
690  this%omega_body = l_body/this%I_body
691  omega = qrot(this%q, this%omega_body)
692 
693  this%vel = this%vel + this%force*dt/(2*this%mass)
694 
695  do i = this%i_start, this%i_stop
696  pos = qrot(this%q, this%pos_body(:,i))
697  ps%vel(:,i) = this%vel + cross(omega, pos)
698  end do
699 
700  end subroutine rigid_body_vv2
701 
702 end module md
subroutine rigid_body_compute_force_torque(this, ps)
Definition: md.f90:592
subroutine, public md_vel(particles, dt)
Definition: md.f90:220
subroutine, public cell_md_pos(cells, particles, dt, md_flag)
Definition: md.f90:80
subroutine rigid_body_vv2(this, ps, dt)
Perform second velocity-Verlet step for quaternion dynamics.
Definition: md.f90:680
subroutine, public cell_md_vel(cells, particles, dt, md_flag)
Definition: md.f90:240
subroutine, public rattle_body_pos(p, links, distances, dt, edges, precision)
Definition: md.f90:322
subroutine rigid_body_init(this, ps, i_start, i_stop, edges)
Definition: md.f90:548
Data for particles.
subroutine, public rattle_dimer_vel(p, d, dt, edges)
Definition: md.f90:299
double precision function, public lj93_zwall(particles, edges, lj_params)
Definition: md.f90:483
pure double precision function, public lj_energy_9_3(z_sq, epsilon, sigma)
pure double precision function, dimension(3), public cross(x1, x2)
Definition: common.f90:405
integer, parameter, public specular_bc
Definition: cell_system.f90:45
subroutine, public rattle_body_vel(p, links, distances, dt, edges, precision)
Definition: md.f90:411
subroutine, public cell_md_pos_ywall(cells, particles, dt, md_flag)
Definition: md.f90:105
subroutine, public apply_pbc(particles, edges)
Definition: md.f90:203
subroutine, public rattle_dimer_pos(p, d, dt, edges)
Definition: md.f90:266
integer, parameter, public wall_bit
Definition: common.f90:51
Utility routines.
Definition: common.f90:12
double precision function, public elastic_network(p, links, distances, k, edges)
Definition: md.f90:517
subroutine, public md_pos(particles, dt)
Definition: md.f90:62
Routines to perform Molecular Dynamics integration.
Definition: md.f90:13
subroutine rigid_body_vv1(this, ps, dt, treshold, edges)
Perform first velocity-Verlet step for quaternion dynamics.
Definition: md.f90:614
pure double precision function, dimension(3), public rel_pos(x, y, L)
Return x-y distance with minimum image convention.
Definition: common.f90:163
pure double precision function, public lj_force_9_3(z, z_sq, epsilon, sigma)
Container for rigid body colloid.
Definition: md.f90:39
Spatial cells.
Definition: cell_system.f90:11
subroutine, public cell_md_pos_zwall(cells, particles, dt, md_flag)
Definition: md.f90:154
Lennard-Jones potential definition.
Definition: interaction.f90:10
integer, parameter, public bounce_back_bc
Definition: cell_system.f90:44