--- trunk/OOPSE/libmdtools/calc_dipole_dipole.F90 2003/04/04 22:22:30 460 +++ trunk/OOPSE/libmdtools/calc_dipole_dipole.F90 2003/04/07 20:06:31 469 @@ -10,10 +10,10 @@ module dipole_dipole #endif implicit none - real(kind=dp), save :: rrf - real(kind=dp), save :: rt - real(kind=dp), save :: rrfsq - real(kind=dp), save :: pre + real(kind=dp), save :: rrf = 0.0 + real(kind=dp), save :: rt = 0.0 + real(kind=dp), save :: rrfsq = 0.0 + real(kind=dp), save :: pre = 0.0 logical, save :: dipole_initialized = .false. contains @@ -39,7 +39,7 @@ contains integer atom1, atom2, me1, me2 real(kind=dp) :: rij, mu1, mu2 - real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5, pre + real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5 real(kind=dp) :: dudx, dudy, dudz, dudu1x, dudu1y, dudu1z real(kind=dp) :: dudu2x, dudu2y, dudu2z, rdotu1, rdotu2, u1dotu2 real(kind=dp) :: taper, dtdr, vterm @@ -81,11 +81,11 @@ contains #endif if( atom1 .eq. 2 )then - write (*,*) 'ul =', ul1(1), ul1(2), ul1(3) + write (0,*) 'ul =', ul1(1), ul1(2), ul1(3) endif if( atom2 .eq. 2 )then - write (*,*) 'ul =', ul2(1), ul2(2), ul2(3) + write (0,*) 'ul =', ul2(1), ul2(2), ul2(3) endif @@ -101,21 +101,20 @@ contains taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3) dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3) endif - + r3 = r2*rij r5 = r3*r2 rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3) rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3) u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3) - + dip2 = pre * mu1 * mu2 - dfact1 = 3.0d0*dip2 / r2 dfact2 = 3.0d0*dip2 / r5 - + vterm = dip2*((u1dotu2/r3) - 3.0d0*(rdotu1*rdotu2/r5)) - + if (do_pot) then #ifdef IS_MPI pot_row(atom1) = pot_row(atom1) + 0.5d0*vterm*taper @@ -139,16 +138,16 @@ contains (5.0d0*(rdotu1*rdotu2)/r5)) - & dfact2*(ul1(3)*rdotu2 + ul2(3)*rdotu1))*taper + & vterm*dtdr*d(3)/rij - + dudu1x = (dip2*((ul2(1)/r3) - (3.0d0*d(1)*rdotu2/r5)))*taper dudu1y = (dip2*((ul2(2)/r3) - (3.0d0*d(2)*rdotu2/r5)))*taper dudu1z = (dip2*((ul2(3)/r3) - (3.0d0*d(3)*rdotu2/r5)))*taper - + dudu2x = (dip2*((ul1(1)/r3) - (3.0d0*d(1)*rdotu1/r5)))*taper dudu2y = (dip2*((ul1(2)/r3) - (3.0d0*d(2)*rdotu1/r5)))*taper dudu2z = (dip2*((ul1(3)/r3) - (3.0d0*d(3)*rdotu1/r5)))*taper - + #ifdef IS_MPI f_Row(1,atom1) = f_Row(1,atom1) - dudx f_Row(2,atom1) = f_Row(2,atom1) - dudy