hdiff output

r31826/chiro.f90 2017-01-26 16:30:09.641146654 +0000 r31825/chiro.f90 2017-01-26 16:30:09.901150168 +0000
  8:         double precision :: r(3), p(3), rm(3, 3), rmd(3, 3, 3)  ! pos, orient, rot mat, rot mat derivs  8:         double precision :: r(3), p(3), rm(3, 3), rmd(3, 3, 3)  ! pos, orient, rot mat, rot mat derivs
  9:         double precision :: mu_hat0(3), mu_hat(3), mu_hat_deriv(3, 3)  9:         double precision :: mu_hat0(3), mu_hat(3), mu_hat_deriv(3, 3)
 10:         integer :: rind, pind   ! indices of values in the big X array 10:         integer :: rind, pind   ! indices of values in the big X array
 11:         ! could put an LJ sites array in here 11:         ! could put an LJ sites array in here
 12:     end type chiro_molecule 12:     end type chiro_molecule
 13:  13: 
 14:     integer :: nmols, x_size    ! number of molecules, size of x array 14:     integer :: nmols, x_size    ! number of molecules, size of x array
 15:  15: 
 16:     type(chiro_molecule), allocatable, target :: molecules(:) 16:     type(chiro_molecule), allocatable, target :: molecules(:)
 17:  17: 
 18:     double precision :: chiro_sigma, chiro_mu, chiro_gamma, chiro_l, jar_radius, k_r 18:     double precision :: chiro_sigma, chiro_mu, chiro_gamma, chiro_l
 19:     double precision :: mu_p2, sin_gamma, cos_gamma 19:     double precision :: mu_p2, sin_gamma, cos_gamma
 20:  20: 
 21: contains 21: contains
 22:     subroutine update_chiro_molecule(mol, r, p, gtest) 22:     subroutine update_chiro_molecule(mol, r, p, gtest)
 23:         use commons, only: periodic, boxlx 23:         use commons, only: periodic, boxlx
 24:         implicit none 24:         implicit none
 25:         type(chiro_molecule), target, intent(inout) :: mol 25:         type(chiro_molecule), target, intent(inout) :: mol
 26:         double precision, intent(inout) :: r(3), p(3) 26:         double precision, intent(inout) :: r(3), p(3)
 27:         logical, intent(in) :: gtest 27:         logical, intent(in) :: gtest
 28:  28: 
 68:         double precision, intent(out) :: dij_p2, xi(3), xj(3), grad(12) 68:         double precision, intent(out) :: dij_p2, xi(3), xj(3), grad(12)
 69:  69: 
 70:         double precision :: rij(3), rij_modp2   ! |rij|^2 70:         double precision :: rij(3), rij_modp2   ! |rij|^2
 71:         double precision :: r_dot_mui, r_dot_muj, mu_dot_mu 71:         double precision :: r_dot_mui, r_dot_muj, mu_dot_mu
 72:         double precision :: cc 72:         double precision :: cc
 73:         double precision :: lambdai, lambdaj ! position of point of contact as measured along rod from the center 73:         double precision :: lambdai, lambdaj ! position of point of contact as measured along rod from the center
 74:         double precision :: offi, offj ! position of naive point of contact as measured from end of rod 74:         double precision :: offi, offj ! position of naive point of contact as measured from end of rod
 75:         double precision :: xij(3) 75:         double precision :: xij(3)
 76:  76: 
 77:         double precision :: dd2_dri(3), dd2_dpi(3), dd2_dpj(3), dd2_dl 77:         double precision :: dd2_dri(3), dd2_dpi(3), dd2_dpj(3), dd2_dl
 78: ! hk286 
 79:         double precision :: dli_dri(3), dli_drj(3), dlj_dri(3), dlj_drj(3) 
 80:         double precision :: dli_dpi(3), dli_dpj(3), dlj_dpi(3), dlj_dpj(3) 
 81:   
 82:         integer :: k 78:         integer :: k
 83:  79: 
 84:         logical :: inti, intj 80:         logical :: inti, intj
 85:  81: 
 86:         ! assume our points are in the interior 82:         ! assume our points are in the interior
 87:         inti = .true. 83:         inti = .true.
 88:         intj = .true. 84:         intj = .true.
 89:  85: 
 90:         rij(:) = ri - rj 86:         rij(:) = ri - rj
 91:         rij_modp2 = dot_product(rij, rij) 87:         rij_modp2 = dot_product(rij, rij)
 92:         r_dot_mui = dot_product(rij, mu_hati) 88:         r_dot_mui = dot_product(rij, mu_hati)
 93:         r_dot_muj = dot_product(rij, mu_hatj) 89:         r_dot_muj = dot_product(rij, mu_hatj)
 94:         mu_dot_mu = dot_product(mu_hati, mu_hatj) 90:         mu_dot_mu = dot_product(mu_hati, mu_hatj)
 95:         cc = 1.d0 - mu_dot_mu * mu_dot_mu 91:         cc = 1.d0 - mu_dot_mu * mu_dot_mu
 96:         92: 
 97:         ! compute naive values for the lambdas 93:         ! compute naive values for the lambdas
 98:         lambdai = (r_dot_mui - mu_dot_mu * r_dot_muj) / cc 94:         lambdai = (r_dot_mui - mu_dot_mu * r_dot_muj) / cc
 99:         lambdaj = (-r_dot_muj + mu_dot_mu * r_dot_mui) / cc 95:         lambdaj = (-r_dot_muj + mu_dot_mu * r_dot_mui) / cc
100:  96: 
101:         if(cc == 0.d0) then 97:         if(cc == 0.d0) then
102:             ! mui is parallel to muj 98:             ! mui is parallel to muj
103:             inti = .false. 99:             inti = .false.
104:             if(r_dot_mui == 0.d0) then100:             if(r_dot_mui == 0.d0) then
105:                 ! rods are side by side: there is no displacement along the101:                 ! rods are side by side: there is no displacement along the
106:                 ! rods. choose the contact point as the middle102:                 ! rods. choose the contact point as the middle
168:             dd2_dri(:) = 2.d0 * xij(:)164:             dd2_dri(:) = 2.d0 * xij(:)
169:             do k = 1, 3165:             do k = 1, 3
170:                 dd2_dpi(k) = dot_product(mui_deriv(k, :), -2.d0 * lambdai * xij(:))166:                 dd2_dpi(k) = dot_product(mui_deriv(k, :), -2.d0 * lambdai * xij(:))
171:                 dd2_dpj(k) = dot_product(muj_deriv(k, :), 2.d0 * lambdaj * xij(:))167:                 dd2_dpj(k) = dot_product(muj_deriv(k, :), 2.d0 * lambdaj * xij(:))
172:             end do168:             end do
173: 169: 
174:             if(inti) then170:             if(inti) then
175:                 ! compute d(d^2)/d(lambda^i)171:                 ! compute d(d^2)/d(lambda^i)
176:                 dd2_dl = 2.d0 * (lambdai - lambdaj * mu_dot_mu - r_dot_mui)172:                 dd2_dl = 2.d0 * (lambdai - lambdaj * mu_dot_mu - r_dot_mui)
177: 173: 
178:                 if (intj) then174:                 ! if at end point, then d(lambda^i)/d(r^i) = 0. if in interior, add on second term.
179:                    ! if at end point, then d(lambda^i)/d(r^i) = 0. if in interior, add on second term.175:                 dd2_dri(:) = dd2_dri(:) + dd2_dl * (mu_hati - mu_dot_mu * mu_hatj) / cc
180:                    dli_dri = (mu_hati - mu_dot_mu * mu_hatj) / cc 
181:                    do k = 1, 3 
182:                       dli_dpi(k) = (2.d0 * mu_dot_mu * dot_product(mui_deriv(k, :), mu_hatj) * & 
183:                            & (r_dot_mui - mu_dot_mu * r_dot_muj) / (cc * cc) + (dot_product(rij, mui_deriv(k, :)) - & 
184:                            & dot_product(mui_deriv(k, :), mu_hatj) * r_dot_muj) / cc ) 
185:                       dli_dpj(k) = (2.d0 * mu_dot_mu * dot_product(mu_hati, muj_deriv(k, :)) * & 
186:                            & (r_dot_mui - mu_dot_mu * r_dot_muj) / (cc * cc) - (dot_product(mu_hati, muj_deriv(k,:)) * & 
187:                            & r_dot_muj + mu_dot_mu * dot_product(rij, muj_deriv(k, :))) / cc ) 
188:                    enddo 
189:                 else 
190:                    dli_dri = mu_hati 
191:                    do k = 1, 3 
192:                       dli_dpi(k) = lambdaj * dot_product(mui_deriv(k, :), mu_hatj) + dot_product(rij, mui_deriv(k, :)) 
193:                       dli_dpj(k) = lambdaj * dot_product(mu_hati, muj_deriv(k, :)) 
194:                    enddo 
195:                 endif 
196:                  
197:                 dd2_dri(:) = dd2_dri(:) + dd2_dl * dli_dri 
198:                 do k = 1, 3176:                 do k = 1, 3
199:                    dd2_dpi(k) = dd2_dpi(k) + dd2_dl * dli_dpi(k)177:                     dd2_dpi(k) = dd2_dpi(k) + dd2_dl * (2.d0 * mu_dot_mu * dot_product(mui_deriv(k, :), mu_hatj) * &
200:                    dd2_dpj(k) = dd2_dpj(k) + dd2_dl * dli_dpj(k)178:                         & (r_dot_mui - mu_dot_mu * r_dot_muj) / (cc * cc) + (dot_product(rij, mui_deriv(k, :)) - &
 179:                         & dot_product(mui_deriv(k, :), mu_hatj) * r_dot_muj) / cc )
 180:                     dd2_dpj(k) = dd2_dpj(k) + dd2_dl * (2.d0 * mu_dot_mu * dot_product(mu_hati, muj_deriv(k, :)) * &
 181:                         & (r_dot_mui - mu_dot_mu * r_dot_muj) / (cc * cc) - (dot_product(mu_hati, muj_deriv(k,:)) * &
 182:                         & r_dot_muj + mu_dot_mu * dot_product(rij, muj_deriv(k, :))) / cc )
201:                 end do183:                 end do
202:              endif184:             end if
203: 185: 
204:             if(intj) then186:             if(intj) then
205:                 ! procedure as above, but for lambda^j187:                 ! procedure as above, but for lambda^j
206:                 dd2_dl = 2.d0 * (lambdaj - lambdai * mu_dot_mu + r_dot_muj)188:                 dd2_dl = 2.d0 * (lambdaj - lambdai * mu_dot_mu + r_dot_muj)
207: 189:                 dd2_dri(:) = dd2_dri(:) + dd2_dl * (-mu_hatj + mu_dot_mu * mu_hati) / cc
208:                 if (inti) then190:                 do k = 1, 3
209:                    dlj_dri = (-mu_hatj + mu_dot_mu * mu_hati) / cc191:                     dd2_dpi(k) = dd2_dpi(k) + dd2_dl * (2.d0 * mu_dot_mu * dot_product(mui_deriv(k, :), mu_hatj) * &
210:                    do k = 1, 3 
211:                       dlj_dpi(k) = (2.d0 * mu_dot_mu * dot_product(mui_deriv(k, :), mu_hatj) * & 
212:                         & (-r_dot_muj + mu_dot_mu * r_dot_mui) / (cc * cc) + (dot_product(mui_deriv(k, :) , mu_hatj) * &192:                         & (-r_dot_muj + mu_dot_mu * r_dot_mui) / (cc * cc) + (dot_product(mui_deriv(k, :) , mu_hatj) * &
213:                         & r_dot_mui + mu_dot_mu * dot_product(rij, mui_deriv(k, :))) / cc )193:                         & r_dot_mui + mu_dot_mu * dot_product(rij, mui_deriv(k, :))) / cc )
214:                       dlj_dpj(k) = (2.d0 * mu_dot_mu * dot_product(mu_hati, muj_deriv(k, :)) * &194:                     dd2_dpj(k) = dd2_dpj(k) + dd2_dl * (2.d0 * mu_dot_mu * dot_product(mu_hati, muj_deriv(k, :)) * &
215:                         & (-r_dot_muj + mu_dot_mu * r_dot_mui) / (cc * cc) + (-dot_product(rij, muj_deriv(k, :)) + &195:                         & (-r_dot_muj + mu_dot_mu * r_dot_mui) / (cc * cc) + (-dot_product(rij, muj_deriv(k, :)) + &
216:                         & dot_product(mu_hati, muj_deriv(k, :)) * r_dot_mui) / cc)196:                         & dot_product(mu_hati, muj_deriv(k, :)) * r_dot_mui) / cc)
217:                    enddo197:                 end do
218:                 else198:             end if
219:                    dlj_dri = - mu_hatj 
220:                    do k = 1, 3 
221:                       dlj_dpi(k) = lambdai * dot_product(mui_deriv(k, :), mu_hatj) 
222:                       dlj_dpj(k) = lambdai * dot_product(mu_hati, muj_deriv(k, :)) - dot_product(rij, muj_deriv(k, :)) 
223:                    enddo 
224:                 endif 
225:                  
226:                 dd2_dri(:) = dd2_dri(:) + dd2_dl * dlj_dri  
227:                 do k = 1, 3 
228:                     dd2_dpi(k) = dd2_dpi(k) + dd2_dl * dlj_dpi(k) 
229:                     dd2_dpj(k) = dd2_dpj(k) + dd2_dl * dlj_dpj(k) 
230:                  end do 
231:               end if 
232:                
233:               grad(1:3) = dd2_dri 
234:               grad(4:6) = -dd2_dri 
235:               grad(7:9) = dd2_dpi 
236:               grad(10:12) = dd2_dpj 
237:            end if 
238: 199: 
239:            return200:             grad(1:3) = dd2_dri
 201:             grad(4:6) = -dd2_dri
 202:             grad(7:9) = dd2_dpi
 203:             grad(10:12) = dd2_dpj
 204:         end if
240: 205: 
 206:         return
241:     end subroutine rod_mindist2207:     end subroutine rod_mindist2
242: 208: 
243:     function max_distance(mols) result(d)209:     function max_distance(mols) result(d)
244:         ! compute the maximum distance of a molecule from the origin210:         ! compute the maximum distance of a molecule from the origin
245: 211: 
246:         use vec3, only: vec_len212:         use vec3, only: vec_len
247: 213: 
248:         implicit none214:         implicit none
249:         type(chiro_molecule), intent(in) :: mols(:)215:         type(chiro_molecule), intent(in) :: mols(:)
250:         double precision :: d216:         double precision :: d
307:     use chiro_module273:     use chiro_module
308:     implicit none274:     implicit none
309:     double precision, intent(inout) :: x(x_size)275:     double precision, intent(inout) :: x(x_size)
310:     logical, intent(in) :: gtest276:     logical, intent(in) :: gtest
311:     double precision, intent(out) :: energy, grad(x_size)277:     double precision, intent(out) :: energy, grad(x_size)
312: 278: 
313:     type(chiro_molecule), pointer :: moli, molj279:     type(chiro_molecule), pointer :: moli, molj
314:     integer moli_ind, molj_ind280:     integer moli_ind, molj_ind
315:     double precision :: grad_pair(12)281:     double precision :: grad_pair(12)
316:     double precision :: energy_contrib, grad_contrib(12)282:     double precision :: energy_contrib, grad_contrib(12)
317:     double precision :: energy_jar, grad_jar(6) 
318: 283: 
319:     double precision :: rij(3), rijmod, rij_hat(3)284:     double precision :: rij(3), rijmod, rij_hat(3)
320:     double precision :: r_pm1, r_pm3, r_pm4, r_pm6, r_pm7, r_pm12, r_pm13 ! r is in reduced units285:     double precision :: r_pm1, r_pm3, r_pm4, r_pm6, r_pm7, r_pm12, r_pm13 ! r is in reduced units
321:     double precision :: mu_dot, mu_cross(3)  ! mu_hati . mu_hatj and mu_hati x mu_hatj286:     double precision :: mu_dot, mu_cross(3)  ! mu_hati . mu_hatj and mu_hati x mu_hatj
322:     integer :: k287:     integer :: k
323: 288: 
324:     double precision :: d_pm1, d_pm2, d_pm6, d_pm12289:     double precision :: d_pm1, d_pm2, d_pm6, d_pm12
325:     double precision :: d2, d2_pm1, d2_pm3, d2_pm6290:     double precision :: d2, d2_pm1, d2_pm3, d2_pm6
326:     double precision :: xi(3), xj(3)291:     double precision :: xi(3), xj(3)
327: 292: 
328:     double precision :: image_moljr(3)293:     double precision :: image_moljr(3)
329:     integer :: lshift, ushift, jshift294:     integer :: lshift, ushift, jshift
330:     logical :: newpair295:     logical :: newpair
331: 296: 
332: ! bjs55 - harmonic force compression for jar potential 
333:  
334:     double precision :: r_xy1, r_xy2 
335:     double precision :: end_point1(3), end_point2(3) 
336:  
337:  
338: ! hk286 
339:     double precision :: epsilon_penalty, delta_penalty 
340:     integer :: power_penalty 
341:  
342:     k_r = 1.0d10 
343:     jar_radius = 5.0 
344:  
345:     epsilon_penalty = 1.0d20 
346:     delta_penalty = 1.0d-3 
347:     power_penalty = 100 
348:  
349:     energy = 0.0297:     energy = 0.0
350:     vt(:) = 0.0298:     vt(:) = 0.0
351:     if(gtest) grad(:) = 0.0299:     if(gtest) grad(:) = 0.0
352: 300: 
353:     ! update the values for all the molecules301:     ! update the values for all the molecules
354:     do moli_ind = 1, nmols302:     do moli_ind = 1, nmols
355:         moli => molecules(moli_ind)303:         moli => molecules(moli_ind)
356:         call update_chiro_molecule(moli, x(moli%rind: moli%rind + 2), x(moli%pind: moli%pind + 2), gtest)304:         call update_chiro_molecule(moli, x(moli%rind: moli%rind + 2), x(moli%pind: moli%pind + 2), gtest)
357:     end do305:     end do
358: 306: 
451: 399: 
452:                 ! U = 4*(d^-12 - d^-6) = 4*((d^2)^-6 - (d^2)^-3)400:                 ! U = 4*(d^-12 - d^-6) = 4*((d^2)^-6 - (d^2)^-3)
453:                 energy_contrib = energy_contrib + 4.d0 * (d2_pm6 - d2_pm3)401:                 energy_contrib = energy_contrib + 4.d0 * (d2_pm6 - d2_pm3)
454:                 if(gtest) then402:                 if(gtest) then
455:                     ! dU/dr^i = 4(-6(d^2)^-7 + 3(d^2)^-4) . d(d^2)/dr^i403:                     ! dU/dr^i = 4(-6(d^2)^-7 + 3(d^2)^-4) . d(d^2)/dr^i
456:                     grad_pair(:) = grad_pair(:) * 4.d0 * (-6.d0 * d2_pm6 + 3.d0 * d2_pm3) * d2_pm1404:                     grad_pair(:) = grad_pair(:) * 4.d0 * (-6.d0 * d2_pm6 + 3.d0 * d2_pm3) * d2_pm1
457:                     grad_contrib = grad_contrib + grad_pair405:                     grad_contrib = grad_contrib + grad_pair
458:                 end if406:                 end if
459: 407: 
460:             end if408:             end if
461:              
462:             ! --- hk286 ------------------------ 
463:             ! --- penalty term contributions --- 
464:             ! compute angular parts only once per set of images 
465:                 mu_dot = dot_product(moli%mu_hat, molj%mu_hat) 
466:                 mu_cross(:) = vec_cross(moli%mu_hat, molj%mu_hat) 
467:  
468:                 rij_hat(:) = rij(:) / rijmod 
469: 409: 
470:                 if ( 1 - mu_dot**2 <  delta_penalty ) then410:             ! add energy contribuions
471:                    energy_contrib = energy_contrib & 
472:                         & + epsilon_penalty * ((1 - mu_dot**2)/delta_penalty-1)**power_penalty  
473:  
474:                    if(gtest) then 
475:                       grad_pair(:) = 0.0 
476:                       do k = 1, 3 
477:                          grad_pair(6 + k) = -2.0d0*power_penalty*epsilon_penalty*mu_dot* & 
478:                               & ((1 - mu_dot**2)/delta_penalty-1)**(power_penalty-1) / delta_penalty * & 
479:                               & dot_product(moli%mu_hat_deriv(k, :), molj%mu_hat) 
480:                           
481:                          grad_pair(9 + k) = -2.0d0*power_penalty*epsilon_penalty*mu_dot* & 
482:                               & ((1 - mu_dot**2)/delta_penalty-1)**(power_penalty-1) / delta_penalty * & 
483:                               & dot_product(moli%mu_hat, molj%mu_hat_deriv(k, :)) 
484:                       end do 
485:                       grad_contrib(:) = grad_contrib(:) + grad_pair(:) 
486:                    end if 
487:                 endif 
488:  
489: ! hk286 
490:  
491:             ! add energy contributions 
492:             energy = energy + energy_contrib411:             energy = energy + energy_contrib
493:             vt(moli_ind) = vt(moli_ind) + energy_contrib412:             vt(moli_ind) = vt(moli_ind) + energy_contrib
494:             vt(molj_ind) = vt(molj_ind) + energy_contrib413:             vt(molj_ind) = vt(molj_ind) + energy_contrib
495: 414: 
496:             ! add gradient contributions415:             ! add gradient contributions
497:             if(gtest) then416:             if(gtest) then
498:                 grad(moli%rind: moli%rind + 2) = grad(moli%rind: moli%rind + 2) + grad_contrib(1: 3)417:                 grad(moli%rind: moli%rind + 2) = grad(moli%rind: moli%rind + 2) + grad_contrib(1: 3)
499:                 grad(molj%rind: molj%rind + 2) = grad(molj%rind: molj%rind + 2) + grad_contrib(4: 6)418:                 grad(molj%rind: molj%rind + 2) = grad(molj%rind: molj%rind + 2) + grad_contrib(4: 6)
500:                 grad(moli%pind: moli%pind + 2) = grad(moli%pind: moli%pind + 2) + grad_contrib(7: 9)419:                 grad(moli%pind: moli%pind + 2) = grad(moli%pind: moli%pind + 2) + grad_contrib(7: 9)
501:                 grad(molj%pind: molj%pind + 2) = grad(molj%pind: molj%pind + 2) + grad_contrib(10: 12)420:                 grad(molj%pind: molj%pind + 2) = grad(molj%pind: molj%pind + 2) + grad_contrib(10: 12)
502:             end if421:             end if
503: 422: 
504:             ! increment pbc shift for inner loop423:             ! increment pbc shift for inner loop
505:             end do424:             end do
506:         end do425:         end do
507: 426: 
508:     end do427:     end do
509: 428: 
510: ! bjs55429:     ! 2D: set all z gradients to zero
511: ! jar harmonic potential terms 
512:  
513:     do moli_ind = 1, nmols 
514:     moli => molecules(moli_ind)   
515:        ! --- jar contributions --- 
516:        end_point1(:) = moli%ri(:) + moli%mu_hati(:) * ( 0.5 * chiro_l ) 
517:        end_point2(:) = moli%ri(:) - moli%mu_hati(:) * ( 0.5 * chiro_l ) 
518:        r_xy1 = sqrt(end_point1(1)**2 + end_point1(2)**2) 
519:        r_xy2 = sqrt(end_point2(1)**2 + end_point2(2)**2) 
520:        energy_jar = 0.0 
521:        grad_jar(:) = 0.0  
522:          
523:        r_xy1 = MAX(r_xy1, r_xy2) 
524:        if ( r_xy1**2 > jar_radius**2 .AND. jar_radius /= 0.d0 ) then 
525:           energy_jar = energy_jar + k_r * (r_xy1 - jar_radius)**2 
526:           energy = energy + energy_jar 
527:           vt(moli_ind) = vt(moli_ind) + energy_jar 
528:           if(gtest) then 
529:              grad_jar(1: 2) = 2 * k_r * end_point1(1: 2) * (1 - ((jar_radius)/r_xy1)) 
530:              grad_jar(3) = 0.0 
531:              do k = 1, 3 
532:                  grad_jar(3 + k) = 2 * k_r * (1 - ((jar_radius)/r_xy1)) * & 
533:                      & (( end_point1(1) * moli%mu_hat_deriv(k, 1)) + ( end_point1(2) * moli%mu_hat_deriv(k, 2))) 
534:              end do 
535:              grad(moli%rind: moli%rind + 2) = grad(moli%rind: moli%rind + 2) + grad_jar(1: 3) 
536:              grad(moli%pind: moli%pind + 2) = grad(moli%pind: moli%pind + 2) + grad_jar(4: 6) 
537:           end if 
538:       end if 
539:     end do 
540:  
541:  ! 2D: set all z gradients to zero 
542:     if(twod) then430:     if(twod) then
543:         do moli_ind = 1, nmols431:         do moli_ind = 1, nmols
544:             grad(molecules(moli_ind)%rind + 2) = 0.0432:             grad(molecules(moli_ind)%rind + 2) = 0.0
545:         end do433:         end do
546:     end if434:     end if
547: 435: 
548: end subroutine chiro436: end subroutine chiro
549: 437: 
550: subroutine chiro_output438: subroutine chiro_output
551:     use commons, only: natoms, nsave, mpit, mynode439:     use commons, only: natoms, nsave, mpit, mynode


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0