hdiff output

r31856/py.f90 2017-01-30 14:30:24.145448146 +0000 r31855/py.f90 2017-01-30 14:30:24.377451264 +0000
  1: module py_module  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/py.f90' in revision 31855
  2: ! 
  3: ! 
  4: ! 
  5: ! sf344> this implementation is still very buggy (GMIN code transported into OPTIM), DO NOT USE!!!! Use the MULTISITEPY keyword instead. 
  6: ! 
  7: ! 
  8: ! 
  9: ! 
 10:  
 11:     ! swo24> equation references are to: 
 12:     !  CW = Chakrabarti & Wales, PCCP 11 1970-1976 (2009) 
 13:     !  PY = Paramonov & Yaliraki, J Chem Phys 123 194111 (2005) 
 14:     ! It would be wise to read and understand these papers before trying to work through this code! 
 15:  
 16:     ! General naming conventions: 
 17:     !  x_p2 = x to power 2 = x ** 2 
 18:     !  x_pm2 = x to power minus 2 = x ** (-2) 
 19:     !  x_mod = modulus of vector x = |x| 
 20:     !  dx_dy = derivative of x with respect to y 
 21:     !  mat_inv = inverse of matrix 'mat' 
 22:  
 23:     implicit none 
 24:      
 25:     double precision, allocatable :: vt(:), coords2(:,:) !, site(:,:)  
 26:     logical, allocatable :: frozen(:), excluded(:) 
 27:     double precision :: step(1), ostep(1), astep(1) 
 28:     double precision :: boxlx, boxly, boxlz 
 29:     logical :: mpit, modulart, percolatet, omove(1), tmove(1) 
 30:     integer :: nsave, mynode, modularcurrentn 
 31:     ! nsave = number of saved minima 
 32:     ! qmin[min] = energy of saved minimum 
 33:     ! qminp[min][xyz] = r and p coordinates of saved minima (as per x array) 
 34:     ! ff[min] = first step when minimum was found 
 35:     double precision, allocatable :: qmin(:), qminp(:,:), ff(:) 
 36:  
 37:     type py_ellipsoid 
 38:         type(py_site), pointer :: site  ! pointer to parent site 
 39:         double precision :: r(3)    ! site position (can vary when using PBCs) 
 40:         double precision :: semiaxes(3), long_semiaxis, short_semiaxis  ! semiaxes of ellipsoid; longest semiaxis 
 41:         ! shape matrix in body frame (i.e., diagonal), in molecule frame, in lab frame, and its inverse (CW eq 23) 
 42:         double precision :: shapemat_bodyframe(3, 3), shapemat_molframe(3, 3), shapemat(3, 3), shapemat_inv(3, 3) 
 43:         double precision :: shapemat_det, shapemat_inv_det  ! determinant and inverse det of shape matrix 
 44:         double precision :: shapemat_inv_cofactors(3, 3)    ! cofactor matrix of inverse shape matrix 
 45:         double precision :: mat_dot_r(3)    ! shape matrix dot position (PY eq 5, second term) 
 46:         double precision :: dF_dr(3), dF_dp(3), dr_dp(3)    ! (CW eq 28); (CW eq 32); (CW eq 12 if at molecule center) 
 47:         double precision :: dmat_dp(3, 3, 3)  ! shape matrix with respect to angle axis (DW eq 33) 
 48:         double precision :: strength    ! coefficient to determiante interaction strength 
 49:     end type py_ellipsoid 
 50:  
 51:     type py_site 
 52:         type(py_molecule), pointer :: molecule  ! pointer to parent molecule 
 53:         type(py_ellipsoid) :: rep, att  ! repulsive and attractive ellipsoids 
 54:         double precision :: r(3), r_molframe(3), p(3), rotmat(3, 3) ! lab and molecule frame positions; rotations 
 55:         logical :: ells_same    ! are the two ellipsoids the same shape and orentiation? 
 56:     end type py_site 
 57:  
 58:     type py_molecule 
 59:         type(py_site), allocatable :: sites(:)  ! container for sites inside the molecule 
 60:         double precision :: r(3), p(3)  ! position and orientation 
 61:         double precision :: rotmat(3, 3), rotmat_deriv(3, 3, 3) ! rotation matrix and derivative 
 62:         integer :: r_index, p_index ! indices of r and p inside the big x array 
 63:     end type py_molecule 
 64:  
 65:     type(py_molecule), target, allocatable :: molecules(:)  ! container for all the molecules 
 66:     integer :: total_num_sites  ! total number of sites in the whole cluster 
 67:  
 68:     ! gcut = value of g (CW eq 25) above which the interaction is cut off 
 69:     double precision :: gcut, gcut_pm2, gcut_pm6, gcut_pm7, gcut_pm8, gcut_pm12, gcut_pm14 
 70:  
 71:     logical :: above_cutoff ! signals between subroutines if ellipsoids are too far apart 
 72:  
 73: contains 
 74:     subroutine initialize_py_potential 
 75:         ! compute values that pertain to the whole PY potential 
 76:  
 77:         use key, only: paramonovcutoff, pcutoff, pysignot 
 78:  
 79:         if(.not. allocated(vt)) allocate(vt(size(molecules)))   ! allocate array for pairwise energies 
 80:  
 81:         if(paramonovcutoff) then 
 82:             gcut = (pcutoff + pysignot) / pysignot 
 83:             gcut_pm2 = 1.d0 / (gcut ** 2) 
 84:             gcut_pm6 = gcut_pm2 ** 3 
 85:             gcut_pm7 = gcut_pm6 / gcut 
 86:             gcut_pm8 = gcut_pm6 * gcut_pm2 
 87:             gcut_pm12 = gcut_pm6 ** 2 
 88:             gcut_pm14 = gcut_pm12 * gcut_pm2 
 89:         end if 
 90:  
 91:     end subroutine initialize_py_potential 
 92:  
 93:     subroutine initialize_py_ellipsoid(site, ell) 
 94:         ! set up the py ellipsoid values for the first time 
 95:  
 96:         implicit none 
 97:         type(py_site), target, intent(in) :: site 
 98:         type(py_ellipsoid), intent(inout) :: ell 
 99:  
100:         double precision :: dummy(3, 3) 
101:         integer :: k 
102:  
103:         ell%site => site    ! set up the parent pointer 
104:         ell%long_semiaxis = maxval(ell%semiaxes) 
105:         ell%short_semiaxis = minval(ell%semiaxes) 
106:  
107:         ell%shapemat_bodyframe(:, :) = 0.d0 
108:         ell%shapemat_det = 1.d0 
109:         do k = 1, 3 
110:             ell%shapemat_bodyframe(k, k) = 1.d0 / (ell%semiaxes(k) ** 2) 
111:             ell%shapemat_det = ell%shapemat_det * ell%shapemat_bodyframe(k, k)  
112:         end do 
113:  
114:         ell%shapemat_inv_det = 1.d0 / ell%shapemat_det 
115:         call rmdrvt(site%p, site%rotmat, dummy, dummy, dummy, .false.) 
116:         ell%shapemat_molframe(:, :) = matmul(site%rotmat, matmul(ell%shapemat_bodyframe, transpose(site%rotmat))) 
117:  
118:     end subroutine initialize_py_ellipsoid 
119:  
120:     subroutine pairwise_py(sitei, sitej, grad_contrib, energy_contrib, gtest) 
121:         ! compute energy and gradient due to a pair of sites 
122:  
123:         use key, only: pysignot, paramonovcutoff, pyepsnot 
124:  
125:         implicit none 
126:         type(py_site), target, intent(inout) :: sitei, sitej 
127:         logical, intent(in) :: gtest 
128:         double precision, intent(out) :: energy_contrib, grad_contrib(12) 
129:  
130:         type(py_ellipsoid), pointer :: elli, ellj   ! pointers to ellipsoids inside the sites 
131:         double precision :: g, g_pm1, g_pm2, g_pm6  ! powers of g 
132:         double precision :: F_pm1_2, dg_dF, dg_dr(3), dg_drmod  ! F_pm1_2 = F ** (-1/2) 
133:         double precision :: e_rep, g_rep(12), e_att, g_att(12) 
134:         double precision :: g_pm7, g_pm12, g_pm13, g_pm14 
135:         double precision :: dU_dF, dU_drmod 
136:  
137:         ! repulsive part 
138:         elli => sitei%rep 
139:         ellj => sitej%rep 
140:         call compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod) 
141:         if(paramonovcutoff .and. above_cutoff) then 
142:             e_rep = 0.d0 
143:             g_rep(:) = 0.d0 
144:         else ! not over the cutoff 
145:             g_pm12 = g_pm6 * g_pm6 
146:             g_pm13 = g_pm12 * g_pm1 
147:  
148:             e_rep = g_pm12 
149:             if(paramonovcutoff) e_rep = e_rep + gcut_pm12 * (6.d0 * gcut_pm2 / g_pm2 - 7.d0) 
150:  
151:             if(gtest) then 
152:                 dU_dF = -2.d0 * g_pm13 * dg_dF ! (CW eq 30) (need 2.d0 for 0.5d0 in dg_dF) 
153:                 dU_drmod = -2.d0 * g_pm13 * (1.d0 - F_pm1_2) / pysignot ! derivate with respect to scalar r_ij 
154:  
155:                 g_rep(1: 3) = -2.d0 * g_pm13 * dg_dr(:)   ! (CW eq 26, but note the error in the paper's sign!) 
156:                 g_rep(7: 9) = dU_dF * elli%dF_dp(:) + dU_drmod * elli%dr_dp(:) ! (CW eq 29, with multisite term) 
157:                 g_rep(10: 12) = dU_dF * ellj%dF_dp(:) + dU_drmod * ellj%dr_dp(:) 
158:  
159:                 if(paramonovcutoff) then 
160:                     ! add splines to achieve cutoff 
161:                     g_rep(1: 3) = g_rep(1: 3) + 2.d0 * gcut_pm14 * g * dg_dr(:) 
162:                     g_rep(7: 9) = g_rep(7: 9) + 2.d0 * gcut_pm14 * g * (dg_dF * elli%dF_dp(:) + dg_drmod * elli%dr_dp(:)) 
163:                     g_rep(10: 12) = g_rep(10: 12) + 2.d0 * gcut_pm14 * g * (dg_dF * ellj%dF_dp(:) + dg_drmod * ellj%dr_dp(:)) 
164:                 end if 
165:  
166:                 g_rep(4: 6) = -g_rep(1: 3)  ! (CW eq 34) 
167:             end if 
168:         end if 
169:  
170:         ! attractive part: if ellipsoids are the same, keep the old values and ellipsoid pointers 
171:         if(.not.(sitei%ells_same .and. sitej%ells_same)) then 
172:             elli => sitei%att 
173:             ellj => sitej%att 
174:             call compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod) 
175:         end if 
176:  
177:         if(paramonovcutoff .and. above_cutoff) then 
178:             e_att = 0.d0 
179:             g_att(:) = 0.d0 
180:              
181:         else ! not over cutoff 
182:             g_pm12 = g_pm6 * g_pm6 
183:             g_pm13 = g_pm12 * g_pm1 
184:  
185:             e_att = -g_pm6 
186:             if(paramonovcutoff) e_att = e_att + gcut_pm6 * (-3.d0 * gcut_pm2 / g_pm2 + 4.d0) 
187:  
188:             if(gtest) then 
189:                 g_pm7 = g_pm6 * g_pm1 
190:                 dU_dF = g_pm7 * dg_dF   ! (CW eq 31, but with factors of 24 and 0.5d0) 
191:                 dU_drmod = g_pm7 * (1.d0 - F_pm1_2) / pysignot 
192:  
193:                 g_att(1: 3) = g_pm7 * dg_dr(:) 
194:                 g_att(7: 9) = dU_dF * elli%dF_dp(:) + dU_drmod * elli%dr_dp(:) 
195:                 g_att(10: 12) = dU_dF * ellj%dF_dp(:) + dU_drmod * ellj%dr_dp(:) 
196:  
197:                 if(paramonovcutoff) then 
198:                     g_att(1: 3) = g_att(1: 3) - gcut_pm8 * g * dg_dr(:) 
199:                     g_att(7: 9) = g_att(7: 9) - gcut_pm8 * g * (dg_dF * elli%dF_dp(:) + dg_drmod * elli%dr_dp(:)) 
200:                     g_att(10: 12) = g_att(10: 12) - gcut_pm8 * g * (dg_dF * ellj%dF_dp(:) + dg_drmod * ellj%dr_dp(:)) 
201:                 end if 
202:  
203:                 g_att(4: 6) = -g_att(1: 3) 
204:             end if 
205:         end if 
206:  
207:         ! add on final factors 
208:         energy_contrib = 4.d0 * pyepsnot * (sitei%rep%strength * sitej%rep%strength * e_rep & 
209:             & + sitei%att%strength * sitej%att%strength * e_att) 
210:         if(gtest) grad_contrib(:) = 24.d0 * pyepsnot * (sitei%rep%strength * sitej%rep%strength * g_rep(:) & 
211:             & + sitei%att%strength * sitej%att%strength * g_att(:)) 
212:     end subroutine pairwise_py 
213:  
214:     subroutine compute_py_values(elli, ellj, gtest, g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr, dg_drmod) 
215:         use key, only: pysignot, paramonovcutoff 
216:         implicit none 
217:  
218:         type(py_ellipsoid), target, intent(inout) :: elli, ellj 
219:         logical, intent(in) :: gtest 
220:         double precision, intent(out) :: g, g_pm1, g_pm2, g_pm6, F_pm1_2, dg_dF, dg_dr(3), dg_drmod 
221:  
222:         type(py_site), pointer :: sitei, sitej 
223:         type(py_molecule), pointer :: moli, molj 
224:         double precision :: lambda_c, F, rij(3), rij_p2, rijmod, rij_hat(3) 
225:         double precision :: x_c(3), xc_minus_ri(3), xc_minus_rj(3)  ! PY eq 9; last term in CW eq 28 
226:         double precision :: xlambdac_term1(3, 3), xlambdac_term2(3) ! PY eq 5 with lambda_c from after PY eq 7 
227:         integer :: i 
228:  
229:         sitei => elli%site 
230:         sitej => ellj%site 
231:         moli => sitei%molecule 
232:         molj => sitej%molecule 
233:  
234:         call compute_contact(elli, ellj, lambda_c, F) 
235:  
236:         rij(:) = elli%r(:) - ellj%r(:) 
237:         rij_p2 = dot_product(rij, rij) 
238:         rijmod = sqrt(rij_p2) 
239:         rij_hat(:) = rij(:) / rijmod 
240:         F_pm1_2 = 1.d0 / sqrt(F) 
241:         g = (rijmod * (1.d0 - F_pm1_2) + pysignot) / pysignot    ! (CW eq 25) 
242:  
243:         if(paramonovcutoff .and. (above_cutoff .or. g > gcut)) then 
244:             ! if over cutoff, do nothing 
245:             above_cutoff = .true. 
246:         else 
247:             g_pm1 = 1.d0 / g 
248:             g_pm2 = g_pm1 * g_pm1 
249:             g_pm6 = g_pm2 * g_pm2 * g_pm2 
250:  
251:             if(gtest) then 
252:                 xlambdac_term1(:, :) = inverse(lambda_c * elli%shapemat + (1.d0 - lambda_c) * ellj%shapemat) 
253:                 xlambdac_term2(:) = lambda_c * elli%mat_dot_r + (1.d0 - lambda_c) * ellj%mat_dot_r 
254:                 x_c(:) = matmul(xlambdac_term1, xlambdac_term2) 
255:                 xc_minus_ri(:) = x_c(:) - elli%r(:) 
256:                 xc_minus_rj(:) = x_c(:) - ellj%r(:) 
257:  
258:                 do i = 1, 3 
259:                     elli%dF_dp(i) = lambda_c * (dot_product(xc_minus_ri, matmul(elli%dmat_dp(i, :, :), xc_minus_ri)) & 
260:                         &- 2.d0 * dot_product(xc_minus_ri, matmul(elli%shapemat, & 
261:                         & matmul(moli%rotmat_deriv(i, :, :), sitei%r_molframe(:)))))    ! CW eq 32 
262:                    ellj%dF_dp(i) = (1.d0 - lambda_c) * (dot_product(xc_minus_rj, matmul(ellj%dmat_dp(i, :, :), xc_minus_rj)) & 
263:                         &- 2.d0 * dot_product(xc_minus_rj, matmul(ellj%shapemat, & 
264:                         & matmul(molj%rotmat_deriv(i, :, :), sitej%r_molframe(:)))))    ! CW eq 35 
265:                 end do 
266:                 elli%dF_dr(:) = -2.d0 * lambda_c * matmul(elli%shapemat, xc_minus_ri)    ! PY eq C2 
267:                 ellj%dF_dr(:) = -elli%dF_dr(:)  ! PY eq C3 
268:  
269:                 dg_dF = 0.5d0 * rijmod * F_pm1_2 / (F * pysignot) 
270:                 dg_dr(:) = (1.d0 - F_pm1_2) * rij_hat(:) / pysignot + dg_dF * elli%dF_dr(:) ! CW eq 17 first term 
271:  
272:                 if(paramonovcutoff) dg_drmod = (1.d0 - F_pm1_2) / pysignot 
273:  
274:                 do i = 1, 3 
275:                     elli%dr_dp(i) = dot_product(rij_hat(:), matmul(moli%rotmat_deriv(i, :, :), sitei%r_molframe(:))) 
276:                     ellj%dr_dp(i) = dot_product(-rij_hat(:), matmul(molj%rotmat_deriv(i, :, :), sitej%r_molframe(:))) 
277:                 end do 
278:             end if 
279:         end if 
280:     end subroutine compute_py_values 
281:  
282:     subroutine compute_contact(elli, ellj, lambda_c, F) 
283:         use key, only: paramonovpbcx, paramonovpbcy, paramonovpbcz, & 
284:             & pysignot, paramonovcutoff, pcutoff, pycfthresh, pycoldfusion 
285:  
286:         implicit none 
287:         type(py_ellipsoid), intent(inout) :: elli, ellj 
288:         double precision, intent(out) :: lambda_c, F 
289:  
290:         integer :: image_list_offset, xyz, o1, o2, image_num 
291:         double precision :: pbc_values(3), image_rj(8, 3), rij(3), rij_p2, rijmod 
292:         double precision :: this_lambda_c, d_R, this_d_R, this_F 
293:  
294:         ! make sure images are initially set properly 
295:         above_cutoff = .false. 
296:         elli%r(:) = elli%site%r(:) 
297:         ellj%r(:) = ellj%site%r(:) 
298:         rij(:) = elli%r(:) - ellj%r(:) 
299:         ellj%mat_dot_r(:) = matmul(ellj%shapemat, ellj%r) 
300:  
301:         if(paramonovpbcx .or. paramonovpbcy .or. paramonovpbcz) then 
302:             ! set up the pbc matrix 
303:             pbc_values(:) = 0.d0 
304:             if(paramonovpbcx) pbc_values(1) = boxlx 
305:             if(paramonovpbcy) pbc_values(2) = boxly 
306:             if(paramonovpbcz) pbc_values(3) = boxlz 
307:  
308:             do xyz = 1, 3 
309:                 if(pbc_values(xyz) /= 0.d0 .and. abs(rij(xyz)) > pbc_values(xyz) / 2.d0) & 
310:                     ellj%r(xyz) = ellj%r(xyz) + sign(pbc_values(xyz), rij(xyz)) 
311:             end do 
312:  
313:             ! redo dependent values, but don't change any orientational stuff 
314:             rij(:) = elli%r(:) - ellj%r(:) 
315:             ellj%mat_dot_r(:) = matmul(ellj%shapemat, ellj%r) 
316:  
317:         end if 
318:  
319:         rijmod = sqrt(dot_product(rij, rij)) 
320:  
321:         if(paramonovcutoff .and. rijmod - elli%long_semiaxis - ellj%long_semiaxis > pcutoff) then 
322:             ! above cutoff 
323:             lambda_c = 0.d0 
324:             F = -1.d0 
325:             above_cutoff = .true. 
326:         else 
327:             call polyecf(elli%shapemat_inv, ellj%shapemat_inv, elli%shapemat_inv_cofactors, & 
328:                 & ellj%shapemat_inv_cofactors, elli%shapemat_inv_det, ellj%shapemat_inv_det, &  
329:                 & rij, lambda_c, F) 
330:                 if(F < pycfthresh) pycoldfusion = .true. 
331:         end if 
332:         ! if ECF falls below given value, then flag for cold fusion of the ellipsoids 
333:  
334:     end subroutine compute_contact 
335:  
336:     subroutine update_py_molecule(mol, r, p, gtest) 
337:         ! changes the r and p of a molecule, then updates all the appropriate mol, site, and ell values 
338:  
339:         use key, only: paramonovpbcx, paramonovpbcy, paramonovpbcz  
340:         implicit none 
341:         double precision, parameter :: pi = 3.14159265358979323846 
342:         type(py_molecule), target, intent(inout) :: mol 
343:         type(py_site), pointer :: site 
344:         double precision, intent(inout) :: r(3), p(3)   ! position and orientation vectors 
345:         logical, intent(in) :: gtest    ! if gradients are required 
346:  
347:         integer :: i 
348:         double precision :: pmod    ! modulus of p 
349:  
350:         ! put molecule back in PBC box 
351:         if(paramonovpbcx) r(1) = r(1) - boxlx * nint(r(1) / boxlx) 
352:         if(paramonovpbcy) r(2) = r(2) - boxly * nint(r(2) / boxly) 
353:         if(paramonovpbcz) r(3) = r(3) - boxlz * nint(r(3) / boxlz) 
354:         mol%r(:) = r(:) 
355:  
356:         ! make sure that 0 < |p| < 2*pi 
357:         pmod = sqrt(dot_product(p, p)) 
358:         if(pmod > 2 * pi) p(:) = p(:) / pmod * mod(pmod, 2 * pi) 
359:         mol%p(:) = p(:) 
360:  
361:         ! recompute the rotation matrices and derivatives 
362:         call rmdrvt(mol%p, mol%rotmat, & 
363:             & mol%rotmat_deriv(1, :, :), mol%rotmat_deriv(2, :, :), mol%rotmat_deriv(3, :, :), gtest) 
364:  
365:         ! loop over sites 
366:         do i = 1, size(mol%sites) 
367:             site => mol%sites(i) 
368:  
369:             site%r(:) = mol%r(:) + matmul(mol%rotmat, site%r_molframe(:))   ! reset site position 
370:             call update_py_ellipsoid(site%rep)  ! update rep ellipsoid 
371:             if(.not. site%ells_same) call update_py_ellipsoid(site%att) ! update att ell if it's different 
372:         end do 
373:  
374:     end subroutine update_py_molecule 
375:  
376:     subroutine update_py_ellipsoid(ell) 
377:         ! update values specific to a ellipsoid 
378:  
379:         implicit none 
380:         type(py_ellipsoid), intent(inout) :: ell 
381:  
382:         type(py_site), pointer :: site 
383:         type(py_molecule), pointer :: mol 
384:         integer :: k 
385:         double precision :: dummy(3, 3) 
386:  
387:         ! identify parent site and molecule 
388:         site => ell%site 
389:         mol => site%molecule 
390:  
391:         ell%r(:) = site%r(:)    ! put ellipsoid at site position 
392:         ell%shapemat(:, :) = matmul(mol%rotmat, matmul(ell%shapemat_molframe, transpose(mol%rotmat))) 
393:         ell%shapemat_inv = adjugate(ell%shapemat) / ell%shapemat_det    ! compute inverse 
394:         ell%shapemat_inv_cofactors = transpose(adjugate(ell%shapemat_inv)) ! cofactor matrix is transpose of adjugate matrix 
395:         ell%mat_dot_r(:) = matmul(ell%shapemat, ell%r) 
396:  
397:         do k = 1, 3 
398:             dummy(:, :) = matmul(mol%rotmat_deriv(k, :, :), ell%shapemat_molframe) 
399:             ell%dmat_dp(k, :, :) = matmul(dummy, transpose(mol%rotmat)) + matmul(mol%rotmat, transpose(dummy)) 
400:         end do 
401:  
402:     end subroutine update_py_ellipsoid 
403:  
404:     function py_overlap(elli, ellj) 
405:         ! return true if two ellipsoids overlap 
406:  
407:         use key, only: pyoverlapthresh, pycoldfusion 
408:         use vec3 
409:         type(py_ellipsoid), intent(inout) :: elli, ellj 
410:         logical :: py_overlap 
411:  
412:         double precision :: rij(3), rijmod, lambda_c, F 
413:  
414:         py_overlap = .false. 
415:         pycoldfusion = .false.  
416:         rij(:) = elli%r(:) - ellj%r(:) 
417:         rijmod = vec_len(rij(:)) 
418:  
419:         ! if ellipsoids could overlap for some orientation 
420:         if(rijmod < elli%short_semiaxis + ellj%short_semiaxis) then 
421:             py_overlap = .true. 
422:         elseif(rijmod < elli%long_semiaxis + ellj%long_semiaxis) then 
423:             call compute_contact(elli, ellj, lambda_c, F) 
424:             if(F < pyoverlapthresh) py_overlap = .true. 
425:         end if 
426:  
427:     end function py_overlap 
428:  
429:     subroutine polyecf(a, b, at, bt, deta, detb, x, x1, smax) 
430:     ! compute ECF using polynomial method and Newton's method to locate the root 
431:  
432:         implicit none 
433:  
434:         ! a = inverse shape matrix of one ellipsoid; at = cofactor matrix of a; deta = det of A 
435:         ! b = other matrix 
436:         ! x[] = separation between ellipsoid centers 
437:         ! x1 = x_n in newton's method, and output value for lambda 
438:         ! smax = s(lambda_c) 
439:         double precision, intent(in) :: a(3,3), at(3,3), b(3,3), bt(3,3) 
440:         double precision, intent(in) :: deta, detb, x(3) 
441:         double precision, intent(out) :: x1, smax 
442:  
443:         ! as = a*, etc 
444:         double precision :: as, bs, cs, ds, es 
445:  
446:         ! x_p2 = x(i) * x(j) terms 
447:         double precision :: x_p2(3, 3) 
448:  
449:         ! f1 = f_1, coefficient of numerator of s(lambda) 
450:         ! g0 = g_0, coefficient of denominator of s(l) 
451:         ! h0 = h_0, coefficient of numerator of s'(l) 
452:         double precision :: f1, f2, f3, f4, g0, g1, g2, g3 
453:         double precision :: h0, h1, h2, h3, h4, h5, h6 
454:  
455:         ! x0 = x_(n-1) in newton's method, test values for lambda 
456:         ! xx = placeholder for powers of x0 (aka lambda) 
457:         ! p  = value of polynomial h at lambda=x0 
458:         ! dp = value of derivative of h 
459:         ! n = numerator of s 
460:         ! d = denominator of s 
461:         double precision :: x0, xx, p, dp, n, d 
462:  
463:         ! i, j = indices 
464:         ! converged = if x1-x0 got sufficiently small 
465:         integer :: i, j 
466:         logical :: converged 
467:  
468:         ! compute x ^ 2 
469:         forall(i = 1:3) 
470:             forall(j = i:3) 
471:                 x_p2(i, j) = x(i) * x(j) 
472:             end forall 
473:         end forall 
474:  
475:         ! compute a* = sum_ij x_i x_j at_ij (noting that at is symmetric) 
476:         as =        x_p2(1, 1)*at(1,1) + x_p2(2, 2)*at(2,2) + x_p2(3, 3)*at(3,3) + & 
477:             & 2.d0*(x_p2(1, 2)*at(1,2) + x_p2(1, 3)*at(1,3) + x_p2(2, 3)*at(2,3))  
478:  
479:         ! compute b* = sum_ij x_i x_j bt_ij (noting that bt is symmetric) 
480:         bs =        x_p2(1, 1)*bt(1,1) + x_p2(2, 2)*bt(2,2) + x_p2(3, 3)*bt(3,3) + & 
481:             & 2.d0*(x_p2(1, 2)*bt(1,2) + x_p2(1, 3)*bt(1,3) + x_p2(2, 3)*bt(2,3)) 
482:  
483:         ! compute c* = sum_ij x_i x_j ct_ij (noting that ct is symmetric) 
484:         cs =       x_p2(1, 1)*(a(2,2)*b(3,3) - a(2,3)*b(3,2) + b(2,2)*a(3,3) - b(2,3)*a(3,2)) + & 
485:            &       x_p2(2, 2)*(a(3,3)*b(1,1) - a(3,1)*b(1,3) + b(3,3)*a(1,1) - b(3,1)*a(1,3)) + & 
486:            &       x_p2(3, 3)*(a(1,1)*b(2,2) - a(1,2)*b(2,1) + b(1,1)*a(2,2) - b(1,2)*a(2,1)) + & 
487:            & 2.d0*(x_p2(1, 2)*(a(2,3)*b(3,1) - a(2,1)*b(3,3) + b(2,3)*a(3,1) - b(2,1)*a(3,3)) + & 
488:            &       x_p2(1, 3)*(a(2,1)*b(3,2) - a(2,2)*b(3,1) + b(2,1)*a(3,2) - b(2,2)*a(3,1)) + & 
489:            &       x_p2(2, 3)*(a(3,1)*b(1,2) - a(3,2)*b(1,1) + b(3,1)*a(1,2) - b(3,2)*a(1,1))) 
490:  
491:         ! compute d* = sum_ij a_ij bt_ij (noting symmetry) 
492:         ds =       a(1,1)*bt(1,1) + a(2,2)*bt(2,2) + a(3,3)*bt(3,3) + & 
493:            & 2.d0*(a(1,2)*bt(1,2) + a(1,3)*bt(1,3) + a(2,3)*bt(2,3)) 
494:  
495:         ! compute e* = sum_ij b_ij at_ij (noting symmetry) 
496:         es =        b(1,1)*at(1,1) + b(2,2)*at(2,2) + b(3,3)*at(3,3) + & 
497:             & 2.d0*(b(1,2)*at(1,2) + b(1,3)*at(1,3) + b(2,3)*at(2,3)) 
498:  
499:         ! compute coefficients of f = numerator of s 
500:         f1 = as 
501:         f2 = -3.d0*as + cs 
502:         f3 = 3.d0*as + bs - 2.d0*cs 
503:         f4 = - as - bs + cs 
504:  
505:         ! compute coefficients of g = denominator of s 
506:         g0 = deta 
507:         g1 = -3.d0*deta + es 
508:         g2 = 3.d0*deta + ds - 2.d0*es 
509:         g3 = -deta + detb - ds + es 
510:  
511:         ! compute coefficients of h = numerator of s' 
512:         h0 = f1*g0 
513:         h1 = 2.d0*f2*g0 
514:         h2 = 3.d0*f3*g0 + f2*g1 - f1*g2 
515:         h3 = 4.d0*f4*g0 + 2.d0*f3*g1 - 2.d0*f1*g3 
516:         h4 = 3.d0*f4*g1 + f3*g2 - f2*g3 
517:         h5 = 2.d0*f4*g2 
518:         h6 = f4*g3 
519:  
520:         ! begin iteration loop. guess initial position. 
521:         x0 = 0.51d0 
522:         converged = .false. 
523:         do i = 1, 25 
524:             ! compute polynomial and its derivative: first, add the constant parts. each 
525:             ! step adds the next power onto both, and xx increases to a new power of x0 
526:             ! each time. 
527:             p  = h0 
528:             dp = h1 
529:             xx = x0 
530:             p  = p  + h1*xx 
531:             dp = dp + 2.d0*h2*xx 
532:             xx = xx*x0 
533:             p  = p  + h2*xx 
534:             dp = dp + 3.d0*h3*xx 
535:             xx = xx*x0 
536:             p  = p  + h3*xx 
537:             dp = dp + 4.d0*h4*xx 
538:             xx = xx*x0 
539:             p  = p  + h4*xx 
540:             dp = dp + 5.d0*h5*xx 
541:             xx = xx*x0 
542:             p  = p  + h5*xx 
543:             dp = dp + 6.d0*h6*xx 
544:             xx = xx*x0 
545:             p  = p  + h6*xx 
546:  
547:             ! compute new position. if dp=0, we are at a critical point and need to  
548:             ! do some extra work. 
549:             if(dp .ne. 0.d0) then 
550:                 ! take the normal newton-rapheson step 
551:                 x1 = x0 - p / dp 
552:  
553:                 ! if x1 ended up outside [0,1], then pretend that x0 was the first point 
554:                 ! computed in a bisection method. so if p>0, go right, and if p<0 go left. 
555:                 if (x1 .lt. 0.d0 .or. x1 .gt. 1.d0) then 
556:                     if(p .gt. 0.d0) x1 = (x0 + 1.d0)/2.d0 
557:                     if(p .lt. 0.d0) x1 =  x0/2.d0 
558:  
559:                 ! if the new position is sufficiently close to the old one, we say we found 
560:                 ! the root. 
561:                 elseif (abs(x1 - x0) .lt. 1d-8) then 
562:                     converged = .true. 
563:                     exit ! leave do loop 
564:                 endif 
565:  
566:             ! if dp=0, take a bisection step. 
567:             else 
568:                 if(p .gt. 0.d0) x1 = (x0 + 1.d0)/2.d0 
569:                 if(p .lt. 0.d0) x1 =  x0/2.d0 
570:             endif 
571:  
572:             ! update old position 
573:             x0 = x1 
574:         enddo 
575:  
576:         ! compute numerator and denominator of s(lambda_c) 
577:         d  = g0 
578:         xx = x1 
579:         n  = f1*xx 
580:         d  = d + g1*xx 
581:         xx = xx*x1 
582:         n  = n + f2*xx 
583:         d  = d + g2*xx 
584:         xx = xx*x1 
585:         n  = n + f3*xx 
586:         d  = d + g3*xx 
587:         xx = xx*x1 
588:         n  = n + f4*xx 
589:  
590:         ! compute s(lambda_c), which gets returned as a parameter 
591:         smax = n / d 
592:  
593:         ! if we did not find a root, notify the user. if this happens, 
594:         ! something has gone very wrong! 
595: !        if (.not. converged) print *, 'polyecf> newton did not converge' 
596:  
597:     end subroutine polyecf 
598:  
599:     function max_distance(mols) result(d) 
600:         ! compute the maximum distance of a molecule from the origin 
601:  
602:         use vec3, only: vec_len 
603:  
604:         implicit none 
605:         type(py_molecule), intent(in) :: mols(:) 
606:         double precision :: d 
607:  
608:         double precision :: this_d 
609:         integer :: i 
610:  
611:         d = 0.d0 
612:  
613:         do i = 1, size(mols) 
614:             this_d = vec_len(mols(i)%r) 
615:             if(this_d > d) d = this_d 
616:         end do 
617:  
618:     end function max_distance 
619:  
620:     function inverse(m) result(inv) 
621:         ! compute inverse of 3x3 matrix 
622:  
623:         implicit none 
624:         double precision, intent(in) :: m(:, :) 
625:         double precision :: inv(3, 3) 
626:  
627:         double precision :: adj(3, 3) 
628:  
629:         adj(:, :) = adjugate(m) 
630:         inv = adj / (m(1, 1) * adj(1, 1) + m(1, 2) * adj(2, 1) + m(1, 3) * adj(3, 1)) 
631:  
632:     end function inverse 
633:  
634:     function adjugate(m) result(a) 
635:         ! return the adjugate (aka adjoint) of a 3x3 matrix 
636:  
637:         ! the adjugate matrix has the properties for 3x3 matrices A: 
638:         !  inverse(A) = adj(A) / det(A) 
639:         !  det(A) = sum_i A_ki adj(A)_ik for k = 1, 2, or 3 
640:  
641:         double precision, intent(in) :: m(3, 3) 
642:         double precision :: a(3, 3) 
643:  
644:         a(1, 1) = m(2, 2) * m(3, 3) - m(3, 2) * m(2, 3) 
645:         a(1, 2) = m(1, 3) * m(3, 2) - m(3, 3) * m(1, 2) 
646:         a(1, 3) = m(1, 2) * m(2, 3) - m(2, 2) * m(1, 3) 
647:         a(2, 1) = m(2, 3) * m(3, 1) - m(3, 3) * m(2, 1) 
648:         a(2, 2) = m(1, 1) * m(3, 3) - m(3, 1) * m(1, 3) 
649:         a(2, 3) = m(1, 3) * m(2, 1) - m(2, 3) * m(1, 1) 
650:         a(3, 1) = m(2, 1) * m(3, 2) - m(3, 1) * m(2, 2) 
651:         a(3, 2) = m(1, 2) * m(3, 1) - m(3, 2) * m(1, 1) 
652:         a(3, 3) = m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2) 
653:  
654:     end function adjugate 
655:  
656: end module py_module 
657:  
658: subroutine py_input 
659:     use commons, only: natoms, nrbsites, rbsite 
660:     use key, only: ntsites, site 
661:     use py_module 
662:     implicit none 
663:     integer, parameter :: py_unit = 299 
664:  
665:     type(py_molecule), pointer :: mol 
666:     type(py_site), pointer :: this_site 
667:     integer :: num_mols, npysites, this_npysite, arity, num_these_mols 
668:     integer :: mol_index, mol_type_index, site_index 
669:     double precision :: number_fraction 
670:     character(len=5) :: label1 
671:     character(len=7) :: label2 
672:     character(len=11) :: label3 
673:  
674:  
675:     num_mols = natoms / 2 
676:     total_num_sites = 0 
677:     allocate(molecules(num_mols)) 
678:     allocate(excluded(num_mols),frozen(natoms)) 
679:     excluded(:)=.false. 
680:     frozen(:)=.false. 
681:     call initialize_py_potential 
682:     !--- parse input 
683:     open(unit=py_unit, file='pysites.xyz', status='old') 
684:     read(py_unit, *) npysites 
685:     if(npysites > 0) then 
686:         ! read in identical particles 
687:         arity = 1 
688:         this_npysite = npysites 
689:         num_these_mols = num_mols 
690:  
691:         ! allocate the SITE array to store the site information for use with rigid body routines 
692:         nrbsites = npysites 
693:     elseif(npysites == 0) then 
694:         ! read in an n-ary mixture 
695:         read(py_unit, *) arity 
696:  
697:         ! pretend that each molecule has only one site, since the rigid body routines can't handle n-ary mixtures 
698:         nrbsites = 1 
699:     end if 
700:  
701:     ! allocate the SITE array for other rigid body routines 
702:     allocate(site(nrbsites,3),rbsite(nrbsites,3)) 
703:     ntsites = num_mols * nrbsites 
704:  
705:     total_num_sites = 0 
706:     mol_index = 1 
707:     do mol_type_index = 1, arity 
708:         read(py_unit, *) ! blank line 
709:         if(npysites == 0) then 
710:             read(py_unit, *) number_fraction   ! read in number fraction for n-ary mixture 
711:             read(py_unit, *) this_npysite      ! read number of site here 
712:  
713:             num_these_mols = nint(num_mols * number_fraction)    ! compute number of this type 
714:         end if 
715:  
716:         if(num_these_mols == 0) then 
717:             read(py_unit, *) ! there will be no molecules of this type, so just toss the information 
718:         else 
719:             ! read in all the sites to the first molecule in this set 
720:             allocate(molecules(mol_index)%sites(this_npysite)) 
721:             do site_index = 1, this_npysite 
722:                 this_site => molecules(mol_index)%sites(site_index) 
723:                 read(py_unit, *) label1, this_site%r_molframe(1), this_site%r_molframe(2), this_site%r_molframe(3), & 
724:                     & label2, this_site%rep%strength, this_site%att%strength, & 
725:                     & this_site%rep%semiaxes(1), this_site%rep%semiaxes(2), this_site%rep%semiaxes(3), & 
726:                     & this_site%att%semiaxes(1), this_site%att%semiaxes(2), this_site%att%semiaxes(3), & 
727:                     & label3, this_site%p(1), this_site%p(2), this_site%p(3) 
728:  
729:                 ! if there is only one type of molecule, then store the site position information in SITE too 
730:                 if(arity == 1) site(site_index, :) = this_site%r_molframe(:) 
731:             end do 
732:             mol_index = mol_index + 1 
733:  
734:             ! copy the first molecule into the rest of the molecules of this type 
735:             do mol_index = mol_index, mol_index + (num_these_mols - 2) 
736:                 allocate(molecules(mol_index)%sites(this_npysite)) 
737:                 molecules(mol_index) = molecules(mol_index - 1) 
738:             end do 
739:  
740:             ! update the total site count 
741:             total_num_sites = total_num_sites + num_these_mols * this_npysite 
742:         end if 
743:     end do 
744:     !--- end parsing input 
745:  
746:     ! initialize 
747:     do mol_index = 1, size(molecules) 
748:         mol => molecules(mol_index) 
749:  
750:         ! initialize molecule values 
751:         mol%r_index = 3 * mol_index - 2 
752:         mol%p_index = mol%r_index + 3 * num_mols 
753:  
754:         do site_index = 1, size(mol%sites) 
755:             this_site => mol%sites(site_index) 
756:             this_site%molecule => mol    ! give all sites a pointer to their parent 
757:  
758:             call initialize_py_ellipsoid(this_site, this_site%rep) 
759:  
760:             if(all(this_site%att%semiaxes == this_site%rep%semiaxes) & 
761:               & .and. (this_site%att%strength == this_site%rep%strength)) then 
762:                 this_site%ells_same = .true. 
763:                 this_site%att = this_site%rep 
764:             else 
765:                 this_site%ells_same = .false. 
766:                 call initialize_py_ellipsoid(this_site, this_site%att) 
767:             end if 
768:                  
769:         end do 
770:     end do 
771: ! sf344> compatibility for the OPTIM code from that written for GMIN: 
772:     rbsite(:,:)=site(:,:) 
773:  
774: end subroutine py_input 
775:  
776: subroutine py(x, grad, energy, gtest) 
777:     use commons, only: natoms  
778:     use key, only: efieldt 
779:     use py_module 
780:     implicit none 
781:     logical :: gtest 
782:     double precision, intent(inout) :: x(3 * natoms) 
783:     double precision, intent(out) :: energy, grad(3 * natoms) 
784:  
785:     type(py_molecule), pointer :: moli, molj 
786:     type(py_site), pointer :: sitei, sitej 
787:     integer :: num_mols, i1 
788:     integer :: moli_index, molj_index, sitei_index, sitej_index 
789:     double precision :: energy_contrib, grad_contrib(12), energy_contrib2, grad_contrib2 
790:  
791:     num_mols = size(molecules) 
792:     gtest = .true. 
793:     energy = 0.d0 
794:     vt(:) = 0.d0 
795:     if(gtest) grad(:) = 0.d0 
796:  
797:     ! update the values for all the molecules 
798:     do moli_index = 1, num_mols 
799:         moli => molecules(moli_index) 
800:         call update_py_molecule(moli, x(moli%r_index: moli%r_index + 2), & 
801:             & x(moli%p_index: moli%p_index + 2), gtest) 
802:     end do 
803:  
804:  
805:     ! outer loop over molecules 
806:     do moli_index = 1, num_mols - 1 
807:         moli => molecules(moli_index) 
808:         ! sf344> do not compute pairwise potential for excluded molecules (modular global optimization) 
809:         if(excluded(moli_index)) cycle 
810:         ! inner loop over molecules 
811:         do molj_index = moli_index + 1, num_mols 
812:             molj => molecules(molj_index) 
813:           ! sf344> do not compute pairwise potential for excluded molecules (modular global optimization) 
814:           if(excluded(molj_index)) cycle 
815:  
816:             ! loop over sites in outer molecule 
817:             do sitei_index = 1, size(moli%sites) 
818:                 sitei => moli%sites(sitei_index) 
819:  
820:                 ! loop over sites in inner molecule 
821:                 do sitej_index = 1, size(molj%sites) 
822:                     sitej => molj%sites(sitej_index) 
823:  
824:                     energy_contrib = 0.d0 
825:                     grad_contrib(:) = 0.d0 
826:  
827:                     ! compute the energy and gradient for this pair of sites 
828:                     call pairwise_py(sitei, sitej, grad_contrib, energy_contrib, gtest) 
829:  
830:                     ! add the energy to total and pairwise 
831:                     energy = energy + energy_contrib 
832:                     vt(moli_index) = vt(moli_index) + energy_contrib 
833:                     vt(molj_index) = vt(molj_index) + energy_contrib 
834:  
835:                     if(gtest) then  ! add gradient contributions 
836:                         grad(moli%r_index: moli%r_index + 2) = grad(moli%r_index: moli%r_index + 2) + grad_contrib(1: 3) 
837:                         grad(molj%r_index: molj%r_index + 2) = grad(molj%r_index: molj%r_index + 2) + grad_contrib(4: 6) 
838:                         grad(moli%p_index: moli%p_index + 2) = grad(moli%p_index: moli%p_index + 2) + grad_contrib(7: 9) 
839:                         grad(molj%p_index: molj%p_index + 2) = grad(molj%p_index: molj%p_index + 2) + grad_contrib(10: 12) 
840:                     end if 
841:  
842:                 end do 
843:             end do 
844:         end do 
845:     end do 
846:     do i1 = 1, natoms/2 
847: !sf344> gravity field with constant gradient, and an r^-12 repulsion around z=0 
848:       if(efieldt) then 
849:         energy_contrib2=0.0 
850:         grad_contrib2=0.0 ! this is a scalar, we have only contribution along the z direction 
851:         call gravity(x(3*i1),energy_contrib2,grad_contrib2) 
852:         energy = energy + energy_contrib2 
853:         grad(3*i1)=grad(3*i1)+grad_contrib2 
854:       end if 
855:     end do 
856:  
857:     ! freeze keyword: set gradients to zero 
858:     if(gtest .and. any(frozen)) then 
859:         do moli_index = 1, num_mols 
860:             moli => molecules(moli_index) 
861:             if(frozen(moli_index)) grad(moli%r_index: moli%r_index + 2) = 0.d0 
862:             if(frozen(moli_index + num_mols)) grad(moli%p_index: moli%p_index + 2) = 0.d0 
863:         end do 
864:     end if 
865:  
866: end subroutine py 
867:  
868: subroutine py_output(output_file,qminp1) 
869:     ! writes ellipsoid.xyz for visualization 
870:  
871:     ! nsave = number of saved minima 
872:     ! qmin[min] = energy of saved minimum 
873:     ! qminp[min][xyz] = r and p coordinates of saved minima (as per x array) 
874:     ! ff[min] = first step when minimum was found 
875:     use py_module 
876:     use commons, only: natoms 
877:     implicit none 
878:     integer, intent(in) :: output_file  ! unit number for output 
879:  
880:     integer :: min_num, mol_num, site_num   ! indices for loops 
881:     type(py_molecule), pointer :: mol       ! pointers for loops 
882:     type(py_site), pointer :: site 
883:     type(py_ellipsoid), pointer :: ell 
884:     double precision :: mat_out(3, 3)   ! rotation matrix used in the output 
885:     character(len=25) :: filename, node, istr 
886:     double precision :: qminp1(3*natoms) 
887:     ! open file for writing 
888: !    if(modulart) then 
889: !         write(istr, *) modularcurrentn 
890: !     if(mpit) then 
891: !         write(node, *) mynode + 1 
892: !         filename = 'ellipsoid_'//trim(adjustl(istr))//'.' // trim(adjustl(node)) // '.xyz' 
893: !     else 
894: !         filename = 'ellipsoid_'//trim(adjustl(istr))//'.xyz' 
895: !     end if 
896: !     open(unit=output_file, file=filename, status='unknown') 
897: !    else 
898: !     if(mpit) then 
899: !         write(node, *) mynode + 1 
900: !         filename = 'ellipsoid.' // trim(adjustl(node)) // '.xyz' 
901: !     else 
902: !         filename = 'ellipsoid.xyz' 
903: !     end if 
904: !     open(unit=output_file, file=filename, status='unknown') 
905: !    end if 
906:     ! loop over saved minima 
907: !    do min_num = 1, nsave 
908:         ! put number of ellipsoids and comment line 
909: !        write(output_file, *) total_num_sites 
910: !        write(output_file, '(A,I10,A,F20.10,A,I10)') 'Energy of minimum', min_num, ' =', qmin(min_num), ' first found at step', ff(min_num) 
911:  
912:         ! loop over molecules 
913:         do mol_num = 1, size(molecules) 
914:             mol => molecules(mol_num) 
915:  
916:             ! use the saved r and p values to recreate the saved configuration 
917:             call update_py_molecule(mol, qminp1(mol%r_index: mol%r_index + 2), & 
918:                 & qminp1(mol%p_index: mol%p_index + 2), .false.) 
919:  
920:             ! loop over sites 
921:             do site_num = 1, size(mol%sites) 
922:                 site => mol%sites(site_num) 
923:                 ell => site%rep     ! use the repulsive ellipsoid for visualization 
924:  
925:                 mat_out = matmul(mol%rotmat, site%rotmat) 
926:  
927:                 write(output_file, '(a5, 2x, 3f20.10, 2x, a8, 12f15.8, 2x, a11, 3f15.8)') & 
928:                     & 'O', site%r(1), site%r(2), site%r(3), & 
929:                     & 'ellipse', 2.d0 * ell%semiaxes(1), 2.d0 * ell%semiaxes(2), 2.d0 * ell%semiaxes(3), & 
930:                     & mat_out(1, 1), mat_out(1, 2), mat_out(1, 3), & 
931:                     & mat_out(2, 1), mat_out(2, 2), mat_out(2, 3), & 
932:                     & mat_out(3, 1), mat_out(3, 2), mat_out(3, 3), & 
933:                     & 'atom_vector', mol%p(1), mol%p(2), mol%p(3)  
934:  
935:             end do 
936:         end do 
937: !    end do 
938:  
939: !    close(output_file) 
940:  
941: end subroutine py_output 
942:  
943: subroutine py_takestep(coords1) 
944:     use commons, only: natoms 
945:     use key, only: myunit, radius, perccut 
946:     use py_module 
947:  
948:     use vec3, only: vec_len, vec_random 
949:     use rotations, only: rot_takestep_aa 
950:  
951:     implicit none 
952:     integer :: np 
953:      
954:     double precision, intent(inout) :: coords1(3*natoms) 
955:     double precision, target :: x(size(coords1(:))) 
956:     double precision, pointer :: r(:), p(:), rsave(:) 
957:     double precision :: step_size, min_vt, dice_roll, stepsave 
958:     type(py_molecule), pointer :: moli, molj 
959:     type(py_site), pointer :: sitei, sitej 
960:     integer :: num_mols, moli_index, molj_index, sitei_index, sitej_index, move_fails 
961:     logical :: stop_pulling ! for angular moves 
962:  
963:     np = 1 
964:     radius=1.0D10 
965:     step(:)=1.0D-1 
966:     astep(:)=1.0D-1 
967:     ostep(:)=1.0D-1 
968:     tmove(:)=.true. 
969:     omove(:)=.true. 
970:     num_mols = size(molecules) 
971:     x(:) = coords1(:) 
972:     min_vt = minval(vt) 
973: !   use stepsave instead of step(np), and reset that for each takestep call, useful for adjusting the step to get non-overlapping moves 
974:     stepsave=step(np) 
975:  
976:     ! update ellipsoids before taking any steps 
977:     do moli_index = 1, num_mols 
978:         moli => molecules(moli_index) 
979:         r => x(moli%r_index: moli%r_index + 2) 
980:         p => x(moli%p_index: moli%p_index + 2) 
981: !        ! sf344> do not compute pairwise potential for excluded molecules (modular global optimization) 
982: !        if(excluded(moli_index)) cycle 
983:  
984:         if(.not. percolatet) then 
985:             if(vec_len(moli%r) > radius) then 
986:                 write(myunit, *) 'py_takestep> initial coord outside container, brining in' 
987:                 r(:) = r(:) - sqrt(radius) * nint(r(:) / sqrt(radius)) 
988:             end if 
989:         end if 
990:  
991:         call update_py_molecule(moli, r, p, .false.) 
992:     end do 
993:  
994:     do moli_index = 1, num_mols 
995:         moli => molecules(moli_index) 
996:         r => x(moli%r_index: moli%r_index + 2) 
997:         p => x(moli%p_index: moli%p_index + 2) 
998:         ! sf344> do not move excluded or frozen molecules (modular global optimization) 
999:         if(excluded(moli_index)) cycle 
1000:         if(frozen(moli_index)) cycle 
1001:  
1002:         ! initialize fail counter. provide landing spot for overlap aborts. if 
1003:         ! the number of fails is too high, just leave this rigid body where it 
1004:         ! is. 
1005:         move_fails = 0 
1006: 95      continue 
1007: !        write(myunit,*) 'move_fails=', move_fails, np, moli_index 
1008: !        write(myunit,*) frozen(moli_index), frozen(molj_index) 
1009:         if(move_fails > 500) STOP   !sf344 debug 
1010: !        if(move_fails > 500) cycle 
1011:  
1012:  
1013:         call update_py_molecule(moli, r, p, .false.) 
1014:  
1015:         ! translational move: uniformly distributed in a sphere of radius step(np) (only if xyz not frozen) 
1016:         if(.not. frozen(moli_index)) then 
1017:             if(tmove(np)) then 
1018:                 if(move_fails>50) stepsave=stepsave*1.1d0 
1019:                 call random_number(step_size) 
1020:                 step_size = sqrt(stepsave * step_size)  ! need to take square root to get uniform volume distribution 
1021:                 r(:) = r(:) + step_size * vec_random() 
1022:             end if 
1023:         end if 
1024:  
1025:         call update_py_molecule(moli, r, p, .false.) 
1026:  
1027:         ! orientational move (only if angle part not frozen) 
1028:         if(.not. frozen(moli_index + num_mols)) then 
1029:             if(omove(np)) then 
1030:                 call random_number(step_size) 
1031:                 call rot_takestep_aa(p(:), step_size * ostep(np)) 
1032:             end if 
1033:         end if 
1034:  
1035:         call update_py_molecule(moli, r, p, .false.) 
1036:         ! angular move: accept step with probability 1 - exp(-beta * delta_energy), where beta = 1/kT = astep 
1037:         ! if xyz or angular parts are frozen, do not move 
1038:         if(.not.(frozen(moli_index) .or. frozen(moli_index + num_mols))) then 
1039:             if(astep(np) > 0.d0) then 
1040:                 call random_number(dice_roll) 
1041:                 if(dice_roll > exp(-astep(np) * (vt(moli_index) - min_vt))) then 
1042:                     r = max_distance(molecules) * vec_random() 
1043:  
1044:                     ! if using percolate, then pull the molecule in until it is within percolation distance of another 
1045:                     if(percolatet) then 
1046:                         stop_pulling = .false. 
1047:                         do 
1048:                             !sf344> bug fix 19.12.2013: molecule coordinates have to be updated, otherwise pulling in has no effect, and this loop can sometimes get infinite! 
1049:                             call update_py_molecule(moli, r, p, .false.) 
1050:                             do molj_index = 1, num_mols 
1051:                                ! sf344> do not move excluded molecules (modular global optimization) 
1052:                                if(excluded(molj_index)) cycle 
1053:                                if(moli_index == molj_index) cycle 
1054:                                 molj => molecules(molj_index) 
1055: !                                    write(myunit, *) 'pulling closer molecules ', moli_index, molj_index, perccut,vec_len(moli%r - molj%r) 
1056:                                 if(vec_len(moli%r - molj%r) < sqrt(perccut)) then  ! perccut = (percolation cutoff)^2! 
1057: !                                    write(myunit, *) 'successfully pulled closer molecules ', moli_index, molj_index, perccut,vec_len(moli%r - molj%r) 
1058:                                     stop_pulling = .true. 
1059:                                     exit 
1060:                                 end if 
1061:                             end do 
1062:  
1063:                             if(stop_pulling) exit 
1064:                             r(:) = 0.95 * r(:) 
1065:                         end do 
1066:                     else ! if using radius, pull in until the molecule is inside the container 
1067:                         do 
1068:                             if(vec_len(r) < radius) exit 
1069:                             r(:) = 0.95 * r(:) 
1070:                         end do 
1071:                     end if 
1072:  
1073:                 end if 
1074:             end if 
1075:         end if 
1076:  
1077:  
1078:         ! update molecule, then check for overlap with other molecules 
1079:         call update_py_molecule(moli, r, p, .false.) 
1080:  
1081:         do molj_index = 1, num_mols 
1082:            ! sf344> do not move excluded or frozen molecules (modular global optimization) 
1083:            if(excluded(molj_index)) cycle 
1084:            if(moli_index == molj_index) cycle 
1085:             molj => molecules(molj_index) 
1086: !            write(myunit,*) 'checking overlap for ', moli_index, molj_index 
1087:  
1088:             do sitei_index = 1, size(moli%sites) 
1089:                 sitei => moli%sites(sitei_index) 
1090:  
1091:                 do sitej_index = 1, size(molj%sites) 
1092:                     sitej => molj%sites(sitej_index) 
1093:  
1094:                     ! if there is overlap, reset the coordinates, increment the 
1095:                     ! fail counter and try again 
1096:                     if(py_overlap(sitei%rep, sitej%rep)) then 
1097: !                        write(myunit,*) 'molecules overlapping: ', moli_index, molj_index 
1098:                         r = coords1(moli%r_index: moli%r_index + 2) 
1099:                         p = coords1(moli%p_index: moli%p_index + 2) 
1100:                         move_fails = move_fails + 1 
1101:                         go to 95 
1102:                     end if 
1103:  
1104:                 end do 
1105:             end do 
1106: !        write(myunit,*) 'distance between molecules ', moli_index, molj_index, vec_len(moli%r - molj%r) 
1107:  
1108:         end do 
1109: !       reset stepsave 
1110:         stepsave=step(np) 
1111: !        print *, 'succeeded with fails: ', move_fails 
1112:     end do 
1113:  
1114:     ! copy local results back 
1115:     coords1(:) = x(:) 
1116:  
1117: end subroutine py_takestep 
1118:  
1119: subroutine py_secder(OLDX,STEST) 
1120: use commons 
1121: !use keyword 
1122: use modhess 
1123: implicit none 
1124:  
1125: DOUBLE PRECISION  :: V(3*NATOMS),EGB,X(3*NATOMS),OLDX(3*NATOMS),VTEMP(2,3*NATOMS),ksi 
1126: INTEGER    :: i,j 
1127: LOGICAL    :: GTEST,STEST 
1128:  
1129:  
1130: ksi=0.00001D0 
1131: X(:)=OLDX(:) 
1132:  
1133: IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS,3*NATOMS)) 
1134: HESS(:,:)=0.0D0 
1135: V(:)=0.0D0 
1136: VTEMP(:,:)=0.0D0 
1137:  
1138: DO i=1,3*NATOMS 
1139:  
1140:          X(i)=X(i)-ksi 
1141:   
1142:          CALL py(X,V,EGB,.true.) 
1143:          VTEMP(1,:)=V(:) 
1144:   
1145:          X(i)=X(i)+2.0D0*ksi 
1146:          CALL py(X,V,EGB,.true.) 
1147:  
1148:          VTEMP(2,:)=V(:) 
1149:  
1150:                 DO j=i,3*NATOMS 
1151:                         HESS(i,j)=(VTEMP(2,j)-VTEMP(1,j))/(2.0D0*ksi) 
1152:                         HESS(j,i)=HESS(i,j) 
1153:                 END DO 
1154: END DO 
1155: end subroutine py_secder 
1156:  
1157: subroutine gravity(x,energy,g) 
1158: use key , only : pygravityc1,pygravityc2 
1159: implicit none 
1160:  
1161: double precision :: x, energy, g 
1162:  
1163: if(x<-9.999D0) then 
1164:         energy = -1.0D10 
1165:         g = -1.0D10 
1166: else 
1167:         energy = pygravityc1/(x+10)**12 + pygravityc2*(x+10) ! using x+10 for backwards compatibility with existing databases  
1168:         g = -12*pygravityc1/(x+10)**13 + pygravityc2    ! so that the potential goes to infinity at z=-10, also shifting the energy range 
1169: end if 
1170: end subroutine gravity 
1171:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0