hdiff output

r31835/chiro.f90 2017-01-27 14:30:34.420277556 +0000 r31834/chiro.f90 2017-01-27 14:30:34.640280529 +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 18:     double precision :: chiro_sigma, chiro_mu, chiro_gamma, chiro_l, jar_radius, k_r
 19:     double precision :: mu_p2, sin_gamma, cos_gamma 19:     double precision :: mu_p2, sin_gamma, cos_gamma
 20:  20: 
 21:  
 22: contains 21: contains
 23:     subroutine update_chiro_molecule(mol, r, p, gtest) 22:     subroutine update_chiro_molecule(mol, r, p, gtest)
 24:         use commons, only: periodic, boxlx 23:         use commons, only: periodic, boxlx
 25:         implicit none 24:         implicit none
 26:         type(chiro_molecule), target, intent(inout) :: mol 25:         type(chiro_molecule), target, intent(inout) :: mol
 27:         double precision, intent(inout) :: r(3), p(3) 26:         double precision, intent(inout) :: r(3), p(3)
 28:         logical, intent(in) :: gtest 27:         logical, intent(in) :: gtest
 29:  28: 
 30:         double precision, parameter :: pi = 3.14159265359 29:         double precision, parameter :: pi = 3.14159265359
 31:         double precision :: pmod 30:         double precision :: pmod
323:     integer :: k322:     integer :: k
324: 323: 
325:     double precision :: d_pm1, d_pm2, d_pm6, d_pm12324:     double precision :: d_pm1, d_pm2, d_pm6, d_pm12
326:     double precision :: d2, d2_pm1, d2_pm3, d2_pm6325:     double precision :: d2, d2_pm1, d2_pm3, d2_pm6
327:     double precision :: xi(3), xj(3)326:     double precision :: xi(3), xj(3)
328: 327: 
329:     double precision :: image_moljr(3)328:     double precision :: image_moljr(3)
330:     integer :: lshift, ushift, jshift329:     integer :: lshift, ushift, jshift
331:     logical :: newpair330:     logical :: newpair
332: 331: 
 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: 
333: ! hk286338: ! hk286
334:     double precision :: epsilon_penalty, delta_penalty339:     double precision :: epsilon_penalty, delta_penalty
335:     integer :: power_penalty340:     integer :: power_penalty
336: 341: 
337:     epsilon_penalty = 1.0d20342:     epsilon_penalty = 1.0d20
338:     delta_penalty = 1.0d-3343:     delta_penalty = 1.0d-3
339:     power_penalty = 100344:     power_penalty = 100
340: 345: 
341:     energy = 0.0346:     energy = 0.0
342:     vt(:) = 0.0347:     vt(:) = 0.0
491:                 grad(molj%rind: molj%rind + 2) = grad(molj%rind: molj%rind + 2) + grad_contrib(4: 6)496:                 grad(molj%rind: molj%rind + 2) = grad(molj%rind: molj%rind + 2) + grad_contrib(4: 6)
492:                 grad(moli%pind: moli%pind + 2) = grad(moli%pind: moli%pind + 2) + grad_contrib(7: 9)497:                 grad(moli%pind: moli%pind + 2) = grad(moli%pind: moli%pind + 2) + grad_contrib(7: 9)
493:                 grad(molj%pind: molj%pind + 2) = grad(molj%pind: molj%pind + 2) + grad_contrib(10: 12)498:                 grad(molj%pind: molj%pind + 2) = grad(molj%pind: molj%pind + 2) + grad_contrib(10: 12)
494:             end if499:             end if
495: 500: 
496:             ! increment pbc shift for inner loop501:             ! increment pbc shift for inner loop
497:             end do502:             end do
498:         end do503:         end do
499: 504: 
500:     end do505:     end do
 506: 
 507: ! bjs55
 508: ! jar harmonic potential terms
 509: 
 510:     do moli_ind = 1, nmols
 511:     moli => molecules(moli_ind)  
 512:        ! --- jar contributions ---
 513:        end_point1(:) = moli%ri(:) + moli%mu_hati(:) * ( 0.5 * chiro_l )
 514:        end_point2(:) = moli%ri(:) - moli%mu_hati(:) * ( 0.5 * chiro_l )
 515:        r_xy1 = sqrt(end_point1(1)**2 + end_point1(2)**2)
 516:        r_xy2 = sqrt(end_point2(1)**2 + end_point2(2)**2)
 517:        energy_jar = 0.0
 518:        grad_jar(:) = 0.0 
 519:         
 520:        r_xy1 = MAX(r_xy1, r_xy2)
 521:        if ( r_xy1**2 > jar_radius**2 .AND. jar_radius /= 0.d0 ) then
 522:           energy_jar = energy_jar + k_r * (r_xy1 - jar_radius)**2
 523:           energy = energy + energy_jar
 524:           vt(moli_ind) = vt(moli_ind) + energy_jar
 525:           if(gtest) then
 526:              grad_jar(1: 2) = 2 * k_r * end_point1(1: 2) * (1 - ((jar_radius)/r_xy1))
 527:              grad_jar(3) = 0.0
 528:              do k = 1, 3
 529:                  grad_jar(3 + k) = 2 * k_r * (1 - ((jar_radius)/r_xy1)) * &
 530:                      & (( end_point1(1) * moli%mu_hat_deriv(k, 1)) + ( end_point1(2) * moli%mu_hat_deriv(k, 2)))
 531:              end do
 532:              grad(moli%rind: moli%rind + 2) = grad(moli%rind: moli%rind + 2) + grad_jar(1: 3)
 533:              grad(moli%pind: moli%pind + 2) = grad(moli%pind: moli%pind + 2) + grad_jar(4: 6)
 534:           end if
 535:       end if
 536:     end do
 537: 
501:  ! 2D: set all z gradients to zero538:  ! 2D: set all z gradients to zero
502:     if(twod) then539:     if(twod) then
503:         do moli_ind = 1, nmols540:         do moli_ind = 1, nmols
504:             grad(molecules(moli_ind)%rind + 2) = 0.0541:             grad(molecules(moli_ind)%rind + 2) = 0.0
505:         end do542:         end do
506:     end if543:     end if
507: 544: 
508: end subroutine chiro545: end subroutine chiro
509: 546: 
510: subroutine chiro_output547: subroutine chiro_output


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0