|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00010 MODULE constraint_4x6 00011 00012 USE atomic_kind_types, ONLY: atomic_kind_type,& 00013 get_atomic_kind 00014 USE constraint_fxd, ONLY: check_fixed_atom_cns_g4x6 00015 USE f77_blas 00016 USE kinds, ONLY: dp 00017 USE linear_systems, ONLY: solve_system 00018 USE mathlib, ONLY: dotprod_3d,& 00019 matvec_3x3 00020 USE molecule_kind_types, ONLY: fixd_constraint_type,& 00021 g4x6_constraint_type,& 00022 get_molecule_kind,& 00023 molecule_kind_type 00024 USE molecule_types_new, ONLY: get_molecule,& 00025 global_constraint_type,& 00026 local_g4x6_constraint_type,& 00027 molecule_type 00028 USE particle_types, ONLY: particle_type 00029 #include "cp_common_uses.h" 00030 00031 IMPLICIT NONE 00032 00033 PRIVATE 00034 PUBLIC :: shake_4x6_int,& 00035 rattle_4x6_int,& 00036 shake_roll_4x6_int,& 00037 rattle_roll_4x6_int,& 00038 shake_4x6_ext,& 00039 rattle_4x6_ext,& 00040 shake_roll_4x6_ext,& 00041 rattle_roll_4x6_ext 00042 00043 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_4x6' 00044 00045 CONTAINS 00046 00047 ! ***************************************************************************** 00052 SUBROUTINE shake_4x6_int( molecule, particle_set, pos, vel, dt, ishake,& 00053 max_sigma, error) 00054 TYPE(molecule_type), POINTER :: molecule 00055 TYPE(particle_type), POINTER :: particle_set( : ) 00056 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00057 REAL(kind=dp), INTENT(in) :: dt 00058 INTEGER, INTENT(IN) :: ishake 00059 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00060 TYPE(cp_error_type), INTENT(inout) :: error 00061 00062 INTEGER :: first_atom, ng4x6 00063 TYPE(fixd_constraint_type), 00064 DIMENSION(:), POINTER :: fixd_list 00065 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00066 TYPE(local_g4x6_constraint_type), 00067 POINTER :: lg4x6( : ) 00068 TYPE(molecule_kind_type), POINTER :: molecule_kind 00069 00070 NULLIFY(fixd_list) 00071 molecule_kind => molecule % molecule_kind 00072 CALL get_molecule_kind ( molecule_kind, ng4x6 = ng4x6,& 00073 fixd_list=fixd_list, g4x6_list = g4x6_list ) 00074 CALL get_molecule ( molecule, first_atom = first_atom, lg4x6=lg4x6 ) 00075 ! Real Shake 00076 CALL shake_4x6_low( fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 00077 particle_set, pos, vel, dt, ishake, max_sigma, error) 00078 00079 END SUBROUTINE shake_4x6_int 00080 00081 ! ***************************************************************************** 00086 SUBROUTINE shake_roll_4x6_int( molecule, particle_set, pos, vel, r_shake, & 00087 dt, ishake, max_sigma, error) 00088 TYPE(molecule_type), POINTER :: molecule 00089 TYPE(particle_type), POINTER :: particle_set( : ) 00090 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00091 REAL(KIND=dp), DIMENSION(:, :), 00092 INTENT(IN) :: r_shake 00093 REAL(kind=dp), INTENT(in) :: dt 00094 INTEGER, INTENT(IN) :: ishake 00095 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00096 TYPE(cp_error_type), INTENT(inout) :: error 00097 00098 INTEGER :: first_atom, ng4x6 00099 TYPE(fixd_constraint_type), 00100 DIMENSION(:), POINTER :: fixd_list 00101 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00102 TYPE(local_g4x6_constraint_type), 00103 POINTER :: lg4x6( : ) 00104 TYPE(molecule_kind_type), POINTER :: molecule_kind 00105 00106 NULLIFY(fixd_list) 00107 molecule_kind => molecule % molecule_kind 00108 CALL get_molecule_kind ( molecule_kind, ng4x6 = ng4x6,& 00109 fixd_list=fixd_list, g4x6_list = g4x6_list ) 00110 CALL get_molecule ( molecule, first_atom = first_atom, lg4x6=lg4x6 ) 00111 ! Real Shake 00112 CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 00113 particle_set, pos, vel, r_shake, dt, ishake, max_sigma,& 00114 error) 00115 00116 END SUBROUTINE shake_roll_4x6_int 00117 00118 ! ***************************************************************************** 00123 SUBROUTINE rattle_4x6_int( molecule, particle_set, vel, dt, error) 00124 TYPE(molecule_type), POINTER :: molecule 00125 TYPE(particle_type), POINTER :: particle_set( : ) 00126 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00127 REAL(kind=dp), INTENT(in) :: dt 00128 TYPE(cp_error_type), INTENT(inout) :: error 00129 00130 INTEGER :: first_atom, ng4x6 00131 TYPE(fixd_constraint_type), 00132 DIMENSION(:), POINTER :: fixd_list 00133 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00134 TYPE(local_g4x6_constraint_type), 00135 POINTER :: lg4x6( : ) 00136 TYPE(molecule_kind_type), POINTER :: molecule_kind 00137 00138 NULLIFY(fixd_list) 00139 molecule_kind => molecule % molecule_kind 00140 CALL get_molecule_kind ( molecule_kind, ng4x6 = ng4x6,& 00141 fixd_list=fixd_list, g4x6_list = g4x6_list ) 00142 CALL get_molecule ( molecule, first_atom = first_atom, lg4x6=lg4x6 ) 00143 ! Real Rattle 00144 CALL rattle_4x6_low( fixd_list, g4x6_list, lg4x6, first_atom, & 00145 particle_set, vel, dt, error) 00146 00147 END SUBROUTINE rattle_4x6_int 00148 00149 ! ***************************************************************************** 00154 SUBROUTINE rattle_roll_4x6_int( molecule, particle_set, vel, r_rattle, dt, veps, error ) 00155 TYPE(molecule_type), POINTER :: molecule 00156 TYPE(particle_type), POINTER :: particle_set( : ) 00157 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00158 REAL(KIND=dp), DIMENSION(:, :), 00159 INTENT(IN) :: r_rattle 00160 REAL(kind=dp), INTENT(in) :: dt 00161 REAL(KIND=dp), DIMENSION(:, :), 00162 INTENT(IN) :: veps 00163 TYPE(cp_error_type), INTENT(inout) :: error 00164 00165 INTEGER :: first_atom, ng4x6 00166 TYPE(fixd_constraint_type), 00167 DIMENSION(:), POINTER :: fixd_list 00168 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00169 TYPE(local_g4x6_constraint_type), 00170 POINTER :: lg4x6( : ) 00171 TYPE(molecule_kind_type), POINTER :: molecule_kind 00172 00173 NULLIFY(fixd_list) 00174 molecule_kind => molecule % molecule_kind 00175 CALL get_molecule_kind ( molecule_kind, ng4x6 = ng4x6,& 00176 fixd_list=fixd_list, g4x6_list = g4x6_list ) 00177 CALL get_molecule ( molecule, first_atom = first_atom, lg4x6=lg4x6 ) 00178 ! Real Rattle 00179 CALL rattle_roll_4x6_low( fixd_list, g4x6_list, lg4x6, first_atom, & 00180 particle_set, vel, r_rattle, dt, veps, error ) 00181 00182 END SUBROUTINE rattle_roll_4x6_int 00183 00184 ! ***************************************************************************** 00189 SUBROUTINE shake_4x6_ext( gci, particle_set, pos, vel, dt, ishake,& 00190 max_sigma, error) 00191 00192 TYPE(global_constraint_type), POINTER :: gci 00193 TYPE(particle_type), POINTER :: particle_set( : ) 00194 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00195 REAL(kind=dp), INTENT(in) :: dt 00196 INTEGER, INTENT(IN) :: ishake 00197 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00198 TYPE(cp_error_type), INTENT(inout) :: error 00199 00200 INTEGER :: first_atom, ng4x6 00201 TYPE(fixd_constraint_type), 00202 DIMENSION(:), POINTER :: fixd_list 00203 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00204 TYPE(local_g4x6_constraint_type), 00205 POINTER :: lg4x6( : ) 00206 00207 first_atom = 1 00208 ng4x6 = gci%ng4x6 00209 fixd_list => gci%fixd_list 00210 g4x6_list => gci%g4x6_list 00211 lg4x6 => gci%lg4x6 00212 ! Real Shake 00213 CALL shake_4x6_low( fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 00214 particle_set, pos, vel, dt, ishake, max_sigma, error) 00215 00216 END SUBROUTINE shake_4x6_ext 00217 00218 ! ***************************************************************************** 00223 SUBROUTINE shake_roll_4x6_ext( gci, particle_set, pos, vel, r_shake, & 00224 dt, ishake, max_sigma, error) 00225 00226 TYPE(global_constraint_type), POINTER :: gci 00227 TYPE(particle_type), POINTER :: particle_set( : ) 00228 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00229 REAL(KIND=dp), DIMENSION(:, :), 00230 INTENT(IN) :: r_shake 00231 REAL(kind=dp), INTENT(in) :: dt 00232 INTEGER, INTENT(IN) :: ishake 00233 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00234 TYPE(cp_error_type), INTENT(inout) :: error 00235 00236 INTEGER :: first_atom, ng4x6 00237 TYPE(fixd_constraint_type), 00238 DIMENSION(:), POINTER :: fixd_list 00239 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00240 TYPE(local_g4x6_constraint_type), 00241 POINTER :: lg4x6( : ) 00242 00243 first_atom = 1 00244 ng4x6 = gci%ng4x6 00245 fixd_list => gci%fixd_list 00246 g4x6_list => gci%g4x6_list 00247 lg4x6 => gci%lg4x6 00248 ! Real Shake 00249 CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 00250 particle_set, pos, vel, r_shake, dt, ishake, max_sigma,& 00251 error) 00252 00253 END SUBROUTINE shake_roll_4x6_ext 00254 00255 ! ***************************************************************************** 00260 SUBROUTINE rattle_4x6_ext( gci, particle_set, vel, dt, error) 00261 00262 TYPE(global_constraint_type), POINTER :: gci 00263 TYPE(particle_type), POINTER :: particle_set( : ) 00264 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00265 REAL(kind=dp), INTENT(in) :: dt 00266 TYPE(cp_error_type), INTENT(inout) :: error 00267 00268 INTEGER :: first_atom, ng4x6 00269 TYPE(fixd_constraint_type), 00270 DIMENSION(:), POINTER :: fixd_list 00271 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00272 TYPE(local_g4x6_constraint_type), 00273 POINTER :: lg4x6( : ) 00274 00275 first_atom = 1 00276 ng4x6 = gci%ng4x6 00277 fixd_list => gci%fixd_list 00278 g4x6_list => gci%g4x6_list 00279 lg4x6 => gci%lg4x6 00280 ! Real Rattle 00281 CALL rattle_4x6_low( fixd_list, g4x6_list, lg4x6, first_atom, & 00282 particle_set, vel, dt, error) 00283 00284 END SUBROUTINE rattle_4x6_ext 00285 00286 ! ***************************************************************************** 00291 SUBROUTINE rattle_roll_4x6_ext( gci, particle_set, vel, r_rattle, dt, veps, error ) 00292 00293 TYPE(global_constraint_type), POINTER :: gci 00294 TYPE(particle_type), POINTER :: particle_set( : ) 00295 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00296 REAL(KIND=dp), DIMENSION(:, :), 00297 INTENT(IN) :: r_rattle 00298 REAL(kind=dp), INTENT(in) :: dt 00299 REAL(KIND=dp), DIMENSION(:, :), 00300 INTENT(IN) :: veps 00301 TYPE(cp_error_type), INTENT(inout) :: error 00302 00303 INTEGER :: first_atom, ng4x6 00304 TYPE(fixd_constraint_type), 00305 DIMENSION(:), POINTER :: fixd_list 00306 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00307 TYPE(local_g4x6_constraint_type), 00308 POINTER :: lg4x6( : ) 00309 00310 first_atom = 1 00311 ng4x6 = gci%ng4x6 00312 fixd_list => gci%fixd_list 00313 g4x6_list => gci%g4x6_list 00314 lg4x6 => gci%lg4x6 00315 ! Real Rattle 00316 CALL rattle_roll_4x6_low( fixd_list, g4x6_list, lg4x6, first_atom, & 00317 particle_set, vel, r_rattle, dt, veps, error ) 00318 00319 END SUBROUTINE rattle_roll_4x6_ext 00320 00321 ! ***************************************************************************** 00325 SUBROUTINE shake_4x6_low( fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 00326 particle_set, pos, vel, dt, ishake, max_sigma, error) 00327 TYPE(fixd_constraint_type), 00328 DIMENSION(:), POINTER :: fixd_list 00329 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00330 TYPE(local_g4x6_constraint_type), 00331 POINTER :: lg4x6( : ) 00332 INTEGER, INTENT(IN) :: ng4x6, first_atom 00333 TYPE(particle_type), POINTER :: particle_set( : ) 00334 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00335 REAL(kind=dp), INTENT(in) :: dt 00336 INTEGER, INTENT(IN) :: ishake 00337 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00338 TYPE(cp_error_type), INTENT(inout) :: error 00339 00340 INTEGER :: iconst, index_a, index_b, 00341 index_c, index_d 00342 REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, 00343 imass1, imass2, imass3, 00344 imass4, sigma 00345 REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, fc4, r0_12, r0_13, r0_14, 00346 r0_23, r0_24, r0_34, r1, r12, r13, r14, r2, r23, r24, r3, r34, r4, v1, 00347 v2, v3, v4, vec 00348 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 00349 REAL(KIND=dp), DIMENSION(6, 6) :: amat, atemp 00350 TYPE(atomic_kind_type), POINTER :: atomic_kind 00351 00352 dtsqby2 = dt*dt*.5_dp 00353 idtsq = 1.0_dp/dt/dt 00354 dtby2 = dt*.5_dp 00355 DO iconst = 1, ng4x6 00356 IF (g4x6_list(iconst)%restraint%active) CYCLE 00357 index_a = g4x6_list ( iconst ) % a + first_atom -1 00358 index_b = g4x6_list ( iconst ) % b + first_atom -1 00359 index_c = g4x6_list ( iconst ) % c + first_atom -1 00360 index_d = g4x6_list ( iconst ) % d + first_atom -1 00361 IF (ishake==1) THEN 00362 r0_12 ( : ) = pos ( :, index_a ) - pos ( :,index_b ) 00363 r0_13 ( : ) = pos ( :, index_a ) - pos ( :,index_c ) 00364 r0_14 ( : ) = pos ( :, index_a ) - pos ( :,index_d ) 00365 r0_23 ( : ) = pos ( :, index_b ) - pos ( :,index_c ) 00366 r0_24 ( : ) = pos ( :, index_b ) - pos ( :,index_d ) 00367 r0_34 ( : ) = pos ( :, index_c ) - pos ( :,index_d ) 00368 atomic_kind => particle_set ( index_a ) % atomic_kind 00369 imass1 = 1.0_dp/atomic_kind%mass 00370 atomic_kind => particle_set ( index_b ) % atomic_kind 00371 imass2 = 1.0_dp/atomic_kind%mass 00372 atomic_kind => particle_set ( index_c ) % atomic_kind 00373 imass3 = 1.0_dp/atomic_kind%mass 00374 atomic_kind => particle_set ( index_d ) % atomic_kind 00375 imass4 = 1.0_dp/atomic_kind%mass 00376 lg4x6 ( iconst ) % fa = -2.0_dp*( lg4x6 ( iconst ) % ra_old - & 00377 lg4x6 ( iconst ) % rb_old ) 00378 lg4x6 ( iconst ) % fb = -2.0_dp*( lg4x6 ( iconst ) % ra_old - & 00379 lg4x6 ( iconst ) % rc_old ) 00380 lg4x6 ( iconst ) % fc = -2.0_dp*( lg4x6 ( iconst ) % ra_old - & 00381 lg4x6 ( iconst ) % rd_old ) 00382 lg4x6 ( iconst ) % fd = -2.0_dp*( lg4x6 ( iconst ) % rb_old - & 00383 lg4x6 ( iconst ) % rc_old ) 00384 lg4x6 ( iconst ) % fe = -2.0_dp*( lg4x6 ( iconst ) % rb_old - & 00385 lg4x6 ( iconst ) % rd_old ) 00386 lg4x6 ( iconst ) % ff = -2.0_dp*( lg4x6 ( iconst ) % rc_old - & 00387 lg4x6 ( iconst ) % rd_old ) 00388 ! Check for fixed atom constraints 00389 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4,& 00390 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 00391 ! construct matrix 00392 amat(1,1)=(imass1+imass2)*DOTPROD_3D(r0_12,lg4x6 ( iconst ) % fa ) 00393 amat(1,2)=imass1*DOTPROD_3D(r0_12, lg4x6 ( iconst) % fb ) 00394 amat(1,3)=imass1*DOTPROD_3D(r0_12, lg4x6 ( iconst ) % fc ) 00395 amat(1,4)=-imass2*DOTPROD_3D(r0_12, lg4x6 ( iconst ) % fd ) 00396 amat(1,5)=-imass2*DOTPROD_3D(r0_12, lg4x6 ( iconst ) % fe ) 00397 amat(1,6)=0.0_dp 00398 00399 amat(2,1)=imass1*DOTPROD_3D(r0_13,lg4x6 ( iconst ) % fa ) 00400 amat(2,2)=(imass1+imass3)*DOTPROD_3D(r0_13, lg4x6 ( iconst) % fb ) 00401 amat(2,3)=imass1*DOTPROD_3D(r0_13, lg4x6 ( iconst ) % fc ) 00402 amat(2,4)=imass3*DOTPROD_3D(r0_13, lg4x6 ( iconst ) % fd ) 00403 amat(2,5)=0.0_dp 00404 amat(2,6)=-imass3*DOTPROD_3D(r0_13, lg4x6 ( iconst ) % ff ) 00405 00406 amat(3,1)=imass1*DOTPROD_3D(r0_14,lg4x6 ( iconst ) % fa ) 00407 amat(3,2)=imass1*DOTPROD_3D(r0_14, lg4x6 ( iconst) % fb ) 00408 amat(3,3)=(imass1+imass4)*DOTPROD_3D(r0_14, lg4x6 ( iconst ) % fc ) 00409 amat(3,4)=0.0_dp 00410 amat(3,5)=imass4*DOTPROD_3D(r0_14, lg4x6 ( iconst ) % fe ) 00411 amat(3,6)=imass4*DOTPROD_3D(r0_14, lg4x6 ( iconst ) % ff ) 00412 00413 amat(4,1)=-imass2*DOTPROD_3D(r0_23,lg4x6 ( iconst ) % fa ) 00414 amat(4,2)=imass3*DOTPROD_3D(r0_23, lg4x6 ( iconst) % fb ) 00415 amat(4,3)=0.0_dp 00416 amat(4,4)=(imass3+imass2)*DOTPROD_3D(r0_23, lg4x6 ( iconst ) % fd ) 00417 amat(4,5)=imass2*DOTPROD_3D(r0_23, lg4x6 ( iconst ) % fe ) 00418 amat(4,6)=-imass3*DOTPROD_3D(r0_23, lg4x6 ( iconst ) % ff ) 00419 00420 amat(5,1)=-imass2*DOTPROD_3D(r0_24,lg4x6 ( iconst ) % fa ) 00421 amat(5,2)=0.0_dp 00422 amat(5,3)=imass4*DOTPROD_3D(r0_24, lg4x6 ( iconst) % fc ) 00423 amat(5,4)=imass2*DOTPROD_3D(r0_24, lg4x6 ( iconst ) % fd ) 00424 amat(5,5)=(imass4+imass2)*DOTPROD_3D(r0_24, lg4x6 ( iconst ) % fe ) 00425 amat(5,6)=imass4*DOTPROD_3D(r0_24, lg4x6 ( iconst ) % ff ) 00426 00427 amat(6,1)=0.0_dp 00428 amat(6,2)=-imass3*DOTPROD_3D(r0_34,lg4x6 ( iconst ) % fb ) 00429 amat(6,3)=imass4*DOTPROD_3D(r0_34, lg4x6 ( iconst) % fc ) 00430 amat(6,4)=-imass3*DOTPROD_3D(r0_34, lg4x6 ( iconst ) % fd ) 00431 amat(6,5)=imass4*DOTPROD_3D(r0_34, lg4x6 ( iconst ) % fe ) 00432 amat(6,6)=(imass3+imass4)*DOTPROD_3D(r0_34, lg4x6 ( iconst ) % ff ) 00433 ! Store values 00434 lg4x6 ( iconst ) % r0_12 = r0_12 00435 lg4x6 ( iconst ) % r0_13 = r0_13 00436 lg4x6 ( iconst ) % r0_14 = r0_14 00437 lg4x6 ( iconst ) % r0_23 = r0_23 00438 lg4x6 ( iconst ) % r0_24 = r0_24 00439 lg4x6 ( iconst ) % r0_34 = r0_34 00440 lg4x6 ( iconst ) % amat = amat 00441 lg4x6 ( iconst ) % imass1 = imass1 00442 lg4x6 ( iconst ) % imass2 = imass2 00443 lg4x6 ( iconst ) % imass3 = imass3 00444 lg4x6 ( iconst ) % imass4 = imass4 00445 lg4x6 ( iconst ) % lambda_old = 0.0_dp 00446 lg4x6 ( iconst ) % del_lambda = 0.0_dp 00447 ELSE 00448 ! Retrieve values 00449 r0_12 = lg4x6 ( iconst ) % r0_12 00450 r0_13 = lg4x6 ( iconst ) % r0_13 00451 r0_14 = lg4x6 ( iconst ) % r0_14 00452 r0_23 = lg4x6 ( iconst ) % r0_23 00453 r0_24 = lg4x6 ( iconst ) % r0_24 00454 r0_34 = lg4x6 ( iconst ) % r0_34 00455 amat = lg4x6 ( iconst ) % amat 00456 imass1 = lg4x6 ( iconst ) % imass1 00457 imass2 = lg4x6 ( iconst ) % imass2 00458 imass3 = lg4x6 ( iconst ) % imass3 00459 imass4 = lg4x6 ( iconst ) % imass4 00460 END IF 00461 00462 ! Iterate until convergence: 00463 vec= lg4x6 ( iconst ) % lambda ( 1 ) * lg4x6 ( iconst ) % fa *(imass1+imass2) + & 00464 lg4x6 ( iconst ) % lambda ( 2 ) * imass1*lg4x6 ( iconst ) % fb + & 00465 lg4x6 ( iconst ) % lambda ( 3 ) * imass1*lg4x6 ( iconst ) % fc - & 00466 lg4x6 ( iconst ) % lambda ( 4 ) * imass2*lg4x6 ( iconst ) % fd - & 00467 lg4x6 ( iconst ) % lambda ( 5 ) * imass2*lg4x6 ( iconst ) % fe 00468 bvec(1,1)=g4x6_list ( iconst ) % dab * g4x6_list ( iconst) % dab & 00469 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_12,r0_12) 00470 vec= lg4x6 ( iconst ) % lambda ( 2 ) * lg4x6 ( iconst ) % fb *(imass1+imass3) + & 00471 lg4x6 ( iconst ) % lambda ( 1 ) * imass1*lg4x6 ( iconst ) % fa + & 00472 lg4x6 ( iconst ) % lambda ( 3 ) * imass1*lg4x6 ( iconst ) % fc + & 00473 lg4x6 ( iconst ) % lambda ( 4 ) * imass3*lg4x6 ( iconst ) % fd - & 00474 lg4x6 ( iconst ) % lambda ( 6 ) * imass3*lg4x6 ( iconst ) % ff 00475 bvec(2,1)=g4x6_list ( iconst ) % dac * g4x6_list ( iconst) % dac & 00476 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_13,r0_13) 00477 vec= lg4x6 ( iconst ) % lambda ( 3 ) * lg4x6 ( iconst ) % fc *(imass1+imass4) + & 00478 lg4x6 ( iconst ) % lambda ( 1 ) * imass1*lg4x6 ( iconst ) % fa + & 00479 lg4x6 ( iconst ) % lambda ( 2 ) * imass1*lg4x6 ( iconst ) % fb + & 00480 lg4x6 ( iconst ) % lambda ( 5 ) * imass4*lg4x6 ( iconst ) % fe + & 00481 lg4x6 ( iconst ) % lambda ( 6 ) * imass4*lg4x6 ( iconst ) % ff 00482 bvec(3,1)=g4x6_list ( iconst ) % dad * g4x6_list ( iconst) % dad & 00483 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_14,r0_14) 00484 vec= lg4x6 ( iconst ) % lambda ( 4 ) * lg4x6 ( iconst ) % fd *(imass2+imass3) - & 00485 lg4x6 ( iconst ) % lambda ( 1 ) * imass2*lg4x6 ( iconst ) % fa + & 00486 lg4x6 ( iconst ) % lambda ( 2 ) * imass3*lg4x6 ( iconst ) % fb + & 00487 lg4x6 ( iconst ) % lambda ( 5 ) * imass2*lg4x6 ( iconst ) % fe - & 00488 lg4x6 ( iconst ) % lambda ( 6 ) * imass3*lg4x6 ( iconst ) % ff 00489 bvec(4,1)=g4x6_list ( iconst ) % dbc * g4x6_list ( iconst) % dbc & 00490 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_23,r0_23) 00491 vec= lg4x6 ( iconst ) % lambda ( 5 ) * lg4x6 ( iconst ) % fe *(imass2+imass4) - & 00492 lg4x6 ( iconst ) % lambda ( 1 ) * imass2*lg4x6 ( iconst ) % fa + & 00493 lg4x6 ( iconst ) % lambda ( 3 ) * imass4*lg4x6 ( iconst ) % fc + & 00494 lg4x6 ( iconst ) % lambda ( 4 ) * imass2*lg4x6 ( iconst ) % fd + & 00495 lg4x6 ( iconst ) % lambda ( 6 ) * imass4*lg4x6 ( iconst ) % ff 00496 bvec(5,1)=g4x6_list ( iconst ) % dbd * g4x6_list ( iconst) % dbd & 00497 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_24,r0_24) 00498 vec= lg4x6 ( iconst ) % lambda ( 6 ) * lg4x6 ( iconst ) % ff *(imass3+imass4) - & 00499 lg4x6 ( iconst ) % lambda ( 2 ) * imass3*lg4x6 ( iconst ) % fb + & 00500 lg4x6 ( iconst ) % lambda ( 3 ) * imass4*lg4x6 ( iconst ) % fc - & 00501 lg4x6 ( iconst ) % lambda ( 4 ) * imass3*lg4x6 ( iconst ) % fd + & 00502 lg4x6 ( iconst ) % lambda ( 5 ) * imass4*lg4x6 ( iconst ) % fe 00503 bvec(6,1)=g4x6_list ( iconst ) % dcd * g4x6_list ( iconst) % dcd & 00504 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_34,r0_34) 00505 bvec = bvec * idtsq 00506 00507 ! get lambda 00508 atemp = amat 00509 CALL solve_system ( atemp, 6, bvec ) 00510 lg4x6 ( iconst ) % lambda ( : ) = bvec ( :, 1 ) 00511 lg4x6 ( iconst ) % del_lambda( : ) = lg4x6 ( iconst ) % lambda ( : ) -& 00512 lg4x6 ( iconst ) % lambda_old ( : ) 00513 lg4x6 ( iconst ) % lambda_old ( : )= lg4x6 ( iconst ) % lambda ( : ) 00514 00515 fc1= lg4x6 ( iconst ) % del_lambda ( 1 ) * lg4x6 ( iconst ) % fa + & 00516 lg4x6 ( iconst ) % del_lambda ( 2 ) * lg4x6 ( iconst ) % fb + & 00517 lg4x6 ( iconst ) % del_lambda ( 3 ) * lg4x6 ( iconst ) % fc 00518 fc2=-lg4x6 ( iconst ) % del_lambda ( 1 ) * lg4x6 ( iconst ) % fa + & 00519 lg4x6 ( iconst ) % del_lambda ( 4 ) * lg4x6 ( iconst ) % fd + & 00520 lg4x6 ( iconst ) % del_lambda ( 5 ) * lg4x6 ( iconst ) % fe 00521 fc3=-lg4x6 ( iconst ) % del_lambda ( 2 ) * lg4x6 ( iconst ) % fb - & 00522 lg4x6 ( iconst ) % del_lambda ( 4 ) * lg4x6 ( iconst ) % fd + & 00523 lg4x6 ( iconst ) % del_lambda ( 6 ) * lg4x6 ( iconst ) % ff 00524 fc4=-lg4x6 ( iconst ) % del_lambda ( 3 ) * lg4x6 ( iconst ) % fc - & 00525 lg4x6 ( iconst ) % del_lambda ( 5 ) * lg4x6 ( iconst ) % fe - & 00526 lg4x6 ( iconst ) % del_lambda ( 6 ) * lg4x6 ( iconst ) % ff 00527 r1 ( : ) = pos ( :, index_a ) + imass1*dtsqby2*fc1 ( : ) 00528 r2 ( : ) = pos ( :, index_b ) + imass2*dtsqby2*fc2 ( : ) 00529 r3 ( : ) = pos ( :, index_c ) + imass3*dtsqby2*fc3 ( : ) 00530 r4 ( : ) = pos ( :, index_d ) + imass4*dtsqby2*fc4 ( : ) 00531 v1 ( : ) = vel ( :, index_a ) + imass1*dtby2*fc1 ( : ) 00532 v2 ( : ) = vel ( :, index_b ) + imass2*dtby2*fc2 ( : ) 00533 v3 ( : ) = vel ( :, index_c ) + imass3*dtby2*fc3 ( : ) 00534 v4 ( : ) = vel ( :, index_d ) + imass4*dtby2*fc4 ( : ) 00535 r12=r1-r2 00536 r13=r1-r3 00537 r14=r1-r4 00538 r23=r2-r3 00539 r24=r2-r4 00540 r34=r3-r4 00541 ! compute the tolerance: 00542 sigma = DOT_PRODUCT(r12,r12) - g4x6_list ( iconst ) % dab * & 00543 g4x6_list ( iconst ) % dab 00544 max_sigma=MAX(max_sigma,ABS(sigma)) 00545 sigma = DOT_PRODUCT(r13,r13) - g4x6_list ( iconst ) % dac * & 00546 g4x6_list ( iconst ) % dac 00547 max_sigma=MAX(max_sigma,ABS(sigma)) 00548 sigma = DOT_PRODUCT(r14,r14) - g4x6_list ( iconst ) % dad * & 00549 g4x6_list ( iconst ) % dad 00550 max_sigma=MAX(max_sigma,ABS(sigma)) 00551 sigma = DOT_PRODUCT(r23,r23) - g4x6_list ( iconst ) % dbc * & 00552 g4x6_list ( iconst ) % dbc 00553 max_sigma=MAX(max_sigma,ABS(sigma)) 00554 sigma = DOT_PRODUCT(r24,r24) - g4x6_list ( iconst ) % dbd * & 00555 g4x6_list ( iconst ) % dbd 00556 max_sigma=MAX(max_sigma,ABS(sigma)) 00557 sigma = DOT_PRODUCT(r34,r34) - g4x6_list ( iconst ) % dcd * & 00558 g4x6_list ( iconst ) % dcd 00559 max_sigma=MAX(max_sigma,ABS(sigma)) 00560 00561 ! update positions with full multiplier 00562 pos ( :, index_a ) = r1 ( : ) 00563 pos ( :, index_b ) = r2 ( : ) 00564 pos ( :, index_c ) = r3 ( : ) 00565 pos ( :, index_d ) = r4 ( : ) 00566 00567 ! update velocites with full multiplier 00568 vel ( :, index_a ) = v1 ( : ) 00569 vel ( :, index_b ) = v2 ( : ) 00570 vel ( :, index_c ) = v3 ( : ) 00571 vel ( :, index_d ) = v4 ( : ) 00572 END DO 00573 END SUBROUTINE shake_4x6_low 00574 00575 ! ***************************************************************************** 00579 SUBROUTINE shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 00580 particle_set, pos, vel, r_shake, dt, ishake, max_sigma,& 00581 error) 00582 TYPE(fixd_constraint_type), 00583 DIMENSION(:), POINTER :: fixd_list 00584 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00585 TYPE(local_g4x6_constraint_type), 00586 POINTER :: lg4x6( : ) 00587 INTEGER, INTENT(IN) :: ng4x6, first_atom 00588 TYPE(particle_type), POINTER :: particle_set( : ) 00589 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00590 REAL(KIND=dp), DIMENSION(:, :), 00591 INTENT(IN) :: r_shake 00592 REAL(kind=dp), INTENT(in) :: dt 00593 INTEGER, INTENT(IN) :: ishake 00594 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00595 TYPE(cp_error_type), INTENT(inout) :: error 00596 00597 INTEGER :: iconst, index_a, index_b, 00598 index_c, index_d 00599 REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, 00600 imass1, imass2, imass3, 00601 imass4, sigma 00602 REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, f_roll4, 00603 f_roll5, f_roll6, fc1, fc2, fc3, fc4, r0_12, r0_13, r0_14, r0_23, 00604 r0_24, r0_34, r1, r12, r13, r14, r2, r23, r24, r3, r34, r4, v1, v2, v3, 00605 v4, vec 00606 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 00607 REAL(KIND=dp), DIMENSION(6, 6) :: amat, atemp 00608 TYPE(atomic_kind_type), POINTER :: atomic_kind 00609 00610 dtsqby2 = dt*dt*.5_dp 00611 idtsq = 1.0_dp/dt/dt 00612 dtby2 = dt*.5_dp 00613 DO iconst = 1, ng4x6 00614 IF (g4x6_list(iconst)%restraint%active) CYCLE 00615 index_a = g4x6_list ( iconst ) % a + first_atom -1 00616 index_b = g4x6_list ( iconst ) % b + first_atom -1 00617 index_c = g4x6_list ( iconst ) % c + first_atom -1 00618 index_d = g4x6_list ( iconst ) % d + first_atom -1 00619 IF (ishake==1) THEN 00620 r0_12 ( : ) = pos ( :, index_a ) - pos ( :,index_b ) 00621 r0_13 ( : ) = pos ( :, index_a ) - pos ( :,index_c ) 00622 r0_23 ( : ) = pos ( :, index_b ) - pos ( :,index_c ) 00623 r0_14 ( : ) = pos ( :, index_a ) - pos ( :,index_d ) 00624 r0_24 ( : ) = pos ( :, index_b ) - pos ( :,index_d ) 00625 r0_34 ( : ) = pos ( :, index_c ) - pos ( :,index_d ) 00626 atomic_kind => particle_set ( index_a ) % atomic_kind 00627 imass1 = 1.0_dp/atomic_kind%mass 00628 atomic_kind => particle_set ( index_b ) % atomic_kind 00629 imass2 = 1.0_dp/atomic_kind%mass 00630 atomic_kind => particle_set ( index_c ) % atomic_kind 00631 imass3 = 1.0_dp/atomic_kind%mass 00632 atomic_kind => particle_set ( index_d ) % atomic_kind 00633 imass4 = 1.0_dp/atomic_kind%mass 00634 lg4x6 ( iconst ) % fa = -2.0_dp*( lg4x6 ( iconst ) % ra_old - & 00635 lg4x6 ( iconst ) % rb_old ) 00636 lg4x6 ( iconst ) % fb = -2.0_dp*( lg4x6 ( iconst ) % ra_old - & 00637 lg4x6 ( iconst ) % rc_old ) 00638 lg4x6 ( iconst ) % fc = -2.0_dp*( lg4x6 ( iconst ) % ra_old - & 00639 lg4x6 ( iconst ) % rd_old ) 00640 lg4x6 ( iconst ) % fd = -2.0_dp*( lg4x6 ( iconst ) % rb_old - & 00641 lg4x6 ( iconst ) % rc_old ) 00642 lg4x6 ( iconst ) % fe = -2.0_dp*( lg4x6 ( iconst ) % rb_old - & 00643 lg4x6 ( iconst ) % rd_old ) 00644 lg4x6 ( iconst ) % ff = -2.0_dp*( lg4x6 ( iconst ) % rc_old - & 00645 lg4x6 ( iconst ) % rd_old ) 00646 ! Check for fixed atom constraints 00647 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4,& 00648 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 00649 ! rotate fconst: 00650 CALL matvec_3x3 (f_roll1, r_shake, lg4x6 ( iconst ) % fa ) 00651 CALL matvec_3x3 (f_roll2, r_shake, lg4x6 ( iconst ) % fb ) 00652 CALL matvec_3x3 (f_roll3, r_shake, lg4x6 ( iconst ) % fc ) 00653 CALL matvec_3x3 (f_roll4, r_shake, lg4x6 ( iconst ) % fd ) 00654 CALL matvec_3x3 (f_roll5, r_shake, lg4x6 ( iconst ) % fe ) 00655 CALL matvec_3x3 (f_roll6, r_shake, lg4x6 ( iconst ) % ff ) 00656 00657 ! construct matrix 00658 amat(1,1)=(imass1+imass2)*DOTPROD_3D(r0_12, f_roll1) 00659 amat(1,2)=imass1*DOTPROD_3D(r0_12, f_roll2 ) 00660 amat(1,3)=imass1*DOTPROD_3D(r0_12, f_roll3 ) 00661 amat(1,4)=-imass2*DOTPROD_3D(r0_12, f_roll4 ) 00662 amat(1,5)=-imass2*DOTPROD_3D(r0_12, f_roll5 ) 00663 amat(1,6)=0.0_dp 00664 00665 amat(2,1)=imass1*DOTPROD_3D(r0_13, f_roll1 ) 00666 amat(2,2)=(imass1+imass3)*DOTPROD_3D(r0_13, f_roll2 ) 00667 amat(2,3)=imass1*DOTPROD_3D(r0_13, f_roll3 ) 00668 amat(2,4)=imass3*DOTPROD_3D(r0_13, f_roll4 ) 00669 amat(2,5)=0.0_dp 00670 amat(2,6)=-imass3*DOTPROD_3D(r0_13, f_roll6 ) 00671 00672 amat(3,1)=imass1*DOTPROD_3D(r0_14, f_roll1 ) 00673 amat(3,2)=imass1*DOTPROD_3D(r0_14, f_roll2 ) 00674 amat(3,3)=(imass1+imass4)*DOTPROD_3D(r0_14, f_roll3 ) 00675 amat(3,4)=0.0_dp 00676 amat(3,5)=imass4*DOTPROD_3D(r0_14, f_roll5 ) 00677 amat(3,6)=imass4*DOTPROD_3D(r0_14, f_roll6 ) 00678 00679 amat(4,1)=-imass2*DOTPROD_3D(r0_23, f_roll1 ) 00680 amat(4,2)=imass3*DOTPROD_3D(r0_23, f_roll2) 00681 amat(4,3)=0.0_dp 00682 amat(4,4)=(imass3+imass2)*DOTPROD_3D(r0_23, f_roll4 ) 00683 amat(4,5)=imass2*DOTPROD_3D(r0_23, f_roll5 ) 00684 amat(4,6)=-imass3*DOTPROD_3D(r0_23, f_roll6 ) 00685 00686 amat(5,1)=-imass2*DOTPROD_3D(r0_24, f_roll1 ) 00687 amat(5,2)=0.0_dp 00688 amat(5,3)=imass4*DOTPROD_3D(r0_24, f_roll3 ) 00689 amat(5,4)=imass2*DOTPROD_3D(r0_24, f_roll4 ) 00690 amat(5,5)=(imass4+imass2)*DOTPROD_3D(r0_24, f_roll5 ) 00691 amat(5,6)=imass4*DOTPROD_3D(r0_24, f_roll6 ) 00692 00693 amat(6,1)=0.0_dp 00694 amat(6,2)=-imass3*DOTPROD_3D(r0_34, f_roll2 ) 00695 amat(6,3)=imass4*DOTPROD_3D(r0_34, f_roll3 ) 00696 amat(6,4)=-imass3*DOTPROD_3D(r0_34, f_roll4 ) 00697 amat(6,5)=imass4*DOTPROD_3D(r0_34, f_roll5 ) 00698 amat(6,6)=(imass3+imass4)*DOTPROD_3D(r0_34, f_roll6 ) 00699 ! Store values 00700 lg4x6 ( iconst ) % r0_12 = r0_12 00701 lg4x6 ( iconst ) % r0_13 = r0_13 00702 lg4x6 ( iconst ) % r0_14 = r0_14 00703 lg4x6 ( iconst ) % r0_23 = r0_23 00704 lg4x6 ( iconst ) % r0_24 = r0_24 00705 lg4x6 ( iconst ) % r0_34 = r0_34 00706 lg4x6 ( iconst ) % f_roll1 = f_roll1 00707 lg4x6 ( iconst ) % f_roll2 = f_roll2 00708 lg4x6 ( iconst ) % f_roll3 = f_roll3 00709 lg4x6 ( iconst ) % f_roll4 = f_roll4 00710 lg4x6 ( iconst ) % f_roll5 = f_roll5 00711 lg4x6 ( iconst ) % f_roll6 = f_roll6 00712 lg4x6 ( iconst ) % amat = amat 00713 lg4x6 ( iconst ) % imass1 = imass1 00714 lg4x6 ( iconst ) % imass2 = imass2 00715 lg4x6 ( iconst ) % imass3 = imass3 00716 lg4x6 ( iconst ) % imass4 = imass4 00717 lg4x6 ( iconst ) % lambda_old = 0.0_dp 00718 lg4x6 ( iconst ) % del_lambda = 0.0_dp 00719 ELSE 00720 ! Retrieve values 00721 r0_12 = lg4x6 ( iconst ) % r0_12 00722 r0_13 = lg4x6 ( iconst ) % r0_13 00723 r0_14 = lg4x6 ( iconst ) % r0_14 00724 r0_23 = lg4x6 ( iconst ) % r0_23 00725 r0_24 = lg4x6 ( iconst ) % r0_24 00726 r0_34 = lg4x6 ( iconst ) % r0_34 00727 f_roll1 = lg4x6 ( iconst ) % f_roll1 00728 f_roll2 = lg4x6 ( iconst ) % f_roll2 00729 f_roll3 = lg4x6 ( iconst ) % f_roll3 00730 f_roll4 = lg4x6 ( iconst ) % f_roll4 00731 f_roll5 = lg4x6 ( iconst ) % f_roll5 00732 f_roll6 = lg4x6 ( iconst ) % f_roll6 00733 amat = lg4x6 ( iconst ) % amat 00734 imass1 = lg4x6 ( iconst ) % imass1 00735 imass2 = lg4x6 ( iconst ) % imass2 00736 imass3 = lg4x6 ( iconst ) % imass3 00737 imass4 = lg4x6 ( iconst ) % imass4 00738 END IF 00739 00740 ! Iterate until convergence: 00741 vec= lg4x6 ( iconst ) % lambda ( 1 ) * f_roll1*(imass1+imass2) + & 00742 lg4x6 ( iconst ) % lambda ( 2 ) * imass1*f_roll2 + & 00743 lg4x6 ( iconst ) % lambda ( 3 ) * imass1*f_roll3 - & 00744 lg4x6 ( iconst ) % lambda ( 4 ) * imass2*f_roll4 - & 00745 lg4x6 ( iconst ) % lambda ( 5 ) * imass2*f_roll5 00746 bvec(1,1)=g4x6_list ( iconst ) % dab * g4x6_list ( iconst) % dab & 00747 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_12,r0_12) 00748 vec= lg4x6 ( iconst ) % lambda ( 2 ) * f_roll2 *(imass1+imass3) + & 00749 lg4x6 ( iconst ) % lambda ( 1 ) * imass1*f_roll1 + & 00750 lg4x6 ( iconst ) % lambda ( 3 ) * imass1*f_roll3 + & 00751 lg4x6 ( iconst ) % lambda ( 4 ) * imass3*f_roll4 - & 00752 lg4x6 ( iconst ) % lambda ( 6 ) * imass3*f_roll6 00753 bvec(2,1)=g4x6_list ( iconst ) % dac * g4x6_list ( iconst) % dac & 00754 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_13,r0_13) 00755 vec= lg4x6 ( iconst ) % lambda ( 3 ) * f_roll3 *(imass1+imass4) + & 00756 lg4x6 ( iconst ) % lambda ( 1 ) * imass1*f_roll1 + & 00757 lg4x6 ( iconst ) % lambda ( 2 ) * imass1*f_roll2 + & 00758 lg4x6 ( iconst ) % lambda ( 5 ) * imass4*f_roll5 + & 00759 lg4x6 ( iconst ) % lambda ( 6 ) * imass4*f_roll6 00760 bvec(3,1)=g4x6_list ( iconst ) % dad * g4x6_list ( iconst) % dad & 00761 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_14,r0_14) 00762 vec= lg4x6 ( iconst ) % lambda ( 4 ) * f_roll4 *(imass2+imass3) - & 00763 lg4x6 ( iconst ) % lambda ( 1 ) * imass2*f_roll1 + & 00764 lg4x6 ( iconst ) % lambda ( 2 ) * imass3*f_roll2 + & 00765 lg4x6 ( iconst ) % lambda ( 5 ) * imass2*f_roll5 - & 00766 lg4x6 ( iconst ) % lambda ( 6 ) * imass3*f_roll6 00767 bvec(4,1)=g4x6_list ( iconst ) % dbc * g4x6_list ( iconst) % dbc & 00768 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_23,r0_23) 00769 vec= lg4x6 ( iconst ) % lambda ( 5 ) * f_roll5 *(imass2+imass4) - & 00770 lg4x6 ( iconst ) % lambda ( 1 ) * imass2*f_roll1 + & 00771 lg4x6 ( iconst ) % lambda ( 3 ) * imass4*f_roll3 + & 00772 lg4x6 ( iconst ) % lambda ( 4 ) * imass2*f_roll4 + & 00773 lg4x6 ( iconst ) % lambda ( 6 ) * imass4*f_roll6 00774 bvec(5,1)=g4x6_list ( iconst ) % dbd * g4x6_list ( iconst) % dbd & 00775 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_24,r0_24) 00776 vec= lg4x6 ( iconst ) % lambda ( 6 ) * f_roll6 *(imass3+imass4) - & 00777 lg4x6 ( iconst ) % lambda ( 2 ) * imass3*f_roll2 + & 00778 lg4x6 ( iconst ) % lambda ( 3 ) * imass4*f_roll3 - & 00779 lg4x6 ( iconst ) % lambda ( 4 ) * imass3*f_roll4 + & 00780 lg4x6 ( iconst ) % lambda ( 5 ) * imass4*f_roll5 00781 bvec(6,1)=g4x6_list ( iconst ) % dcd * g4x6_list ( iconst) % dcd & 00782 -dtsqby2*dtsqby2*DOTPROD_3D(vec,vec)-DOTPROD_3D(r0_34,r0_34) 00783 bvec = bvec * idtsq 00784 00785 ! get lambda 00786 atemp = amat 00787 CALL solve_system ( atemp, 6, bvec ) 00788 lg4x6 ( iconst ) % lambda ( : ) = bvec ( :, 1 ) 00789 lg4x6 ( iconst ) % del_lambda( : ) = lg4x6 ( iconst ) % lambda ( : ) -& 00790 lg4x6 ( iconst ) % lambda_old ( : ) 00791 lg4x6 ( iconst ) % lambda_old ( : )= lg4x6 ( iconst ) % lambda ( : ) 00792 00793 fc1= lg4x6 ( iconst ) % del_lambda ( 1 ) * lg4x6 ( iconst ) % fa + & 00794 lg4x6 ( iconst ) % del_lambda ( 2 ) * lg4x6 ( iconst ) % fb + & 00795 lg4x6 ( iconst ) % del_lambda ( 3 ) * lg4x6 ( iconst ) % fc 00796 fc2=-lg4x6 ( iconst ) % del_lambda ( 1 ) * lg4x6 ( iconst ) % fa + & 00797 lg4x6 ( iconst ) % del_lambda ( 4 ) * lg4x6 ( iconst ) % fd + & 00798 lg4x6 ( iconst ) % del_lambda ( 5 ) * lg4x6 ( iconst ) % fe 00799 fc3=-lg4x6 ( iconst ) % del_lambda ( 2 ) * lg4x6 ( iconst ) % fb - & 00800 lg4x6 ( iconst ) % del_lambda ( 4 ) * lg4x6 ( iconst ) % fd + & 00801 lg4x6 ( iconst ) % del_lambda ( 6 ) * lg4x6 ( iconst ) % ff 00802 fc4=-lg4x6 ( iconst ) % del_lambda ( 3 ) * lg4x6 ( iconst ) % fc - & 00803 lg4x6 ( iconst ) % del_lambda ( 5 ) * lg4x6 ( iconst ) % fe - & 00804 lg4x6 ( iconst ) % del_lambda ( 6 ) * lg4x6 ( iconst ) % ff 00805 CALL MATVEC_3x3(vec,r_shake,fc1) 00806 r1 ( : ) = pos ( :, index_a ) + imass1*dtsqby2*vec 00807 CALL MATVEC_3x3(vec,r_shake,fc2) 00808 r2 ( : ) = pos ( :, index_b ) + imass2*dtsqby2*vec 00809 CALL MATVEC_3x3(vec,r_shake,fc3) 00810 r3 ( : ) = pos ( :, index_c ) + imass3*dtsqby2*vec 00811 CALL MATVEC_3x3(vec,r_shake,fc4) 00812 r4 ( : ) = pos ( :, index_d ) + imass4*dtsqby2*vec 00813 CALL MATVEC_3x3(vec,r_shake,fc1) 00814 v1 ( : ) = vel ( :, index_a ) + imass1*dtby2*vec 00815 CALL MATVEC_3x3(vec,r_shake,fc2) 00816 v2 ( : ) = vel ( :, index_b ) + imass2*dtby2*vec 00817 CALL MATVEC_3x3(vec,r_shake,fc3) 00818 v3 ( : ) = vel ( :, index_c ) + imass3*dtby2*vec 00819 CALL MATVEC_3x3(vec,r_shake,fc4) 00820 v4 ( : ) = vel ( :, index_d ) + imass4*dtby2*vec 00821 00822 r12=r1-r2 00823 r13=r1-r3 00824 r23=r2-r3 00825 r14=r1-r4 00826 r24=r2-r4 00827 r34=r3-r4 00828 ! compute the tolerance: 00829 sigma = DOT_PRODUCT(r12,r12) - g4x6_list ( iconst ) % dab * & 00830 g4x6_list ( iconst ) % dab 00831 max_sigma=MAX(max_sigma,ABS(sigma)) 00832 sigma = DOT_PRODUCT(r13,r13) - g4x6_list ( iconst ) % dac * & 00833 g4x6_list ( iconst ) % dac 00834 max_sigma=MAX(max_sigma,ABS(sigma)) 00835 sigma = DOT_PRODUCT(r14,r14) - g4x6_list ( iconst ) % dad * & 00836 g4x6_list ( iconst ) % dad 00837 max_sigma=MAX(max_sigma,ABS(sigma)) 00838 sigma = DOT_PRODUCT(r23,r23) - g4x6_list ( iconst ) % dbc * & 00839 g4x6_list ( iconst ) % dbc 00840 max_sigma=MAX(max_sigma,ABS(sigma)) 00841 sigma = DOT_PRODUCT(r24,r24) - g4x6_list ( iconst ) % dbd * & 00842 g4x6_list ( iconst ) % dbd 00843 max_sigma=MAX(max_sigma,ABS(sigma)) 00844 sigma = DOT_PRODUCT(r34,r34) - g4x6_list ( iconst ) % dcd * & 00845 g4x6_list ( iconst ) % dcd 00846 max_sigma=MAX(max_sigma,ABS(sigma)) 00847 00848 ! update positions with full multiplier 00849 pos ( :, index_a ) = r1 ( : ) 00850 pos ( :, index_b ) = r2 ( : ) 00851 pos ( :, index_c ) = r3 ( : ) 00852 pos ( :, index_d ) = r4 ( : ) 00853 00854 ! update velocites with full multiplier 00855 vel ( :, index_a ) = v1 ( : ) 00856 vel ( :, index_b ) = v2 ( : ) 00857 vel ( :, index_c ) = v3 ( : ) 00858 vel ( :, index_d ) = v4 ( : ) 00859 END DO 00860 END SUBROUTINE shake_roll_4x6_low 00861 00862 ! ***************************************************************************** 00866 SUBROUTINE rattle_4x6_low( fixd_list, g4x6_list, lg4x6, first_atom, & 00867 particle_set, vel, dt, error) 00868 00869 TYPE(fixd_constraint_type), 00870 DIMENSION(:), POINTER :: fixd_list 00871 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 00872 TYPE(local_g4x6_constraint_type), 00873 POINTER :: lg4x6( : ) 00874 INTEGER, INTENT(IN) :: first_atom 00875 TYPE(particle_type), POINTER :: particle_set( : ) 00876 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00877 REAL(kind=dp), INTENT(in) :: dt 00878 TYPE(cp_error_type), INTENT(inout) :: error 00879 00880 INTEGER :: iconst, index_a, index_b, 00881 index_c, index_d 00882 REAL(KIND=dp) :: dtby2, idt, imass1, imass2, 00883 imass3, imass4, mass 00884 REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, fc4, r12, r13, 00885 r14, r23, r24, r34, v12, v13, 00886 v14, v23, v24, v34 00887 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 00888 REAL(KIND=dp), DIMENSION(6, 6) :: amat 00889 TYPE(atomic_kind_type), POINTER :: atomic_kind 00890 00891 idt = 1.0_dp/dt 00892 dtby2 = dt*.5_dp 00893 DO iconst = 1, SIZE ( g4x6_list ) 00894 IF (g4x6_list(iconst)%restraint%active) CYCLE 00895 index_a = g4x6_list ( iconst ) % a + first_atom -1 00896 index_b = g4x6_list ( iconst ) % b + first_atom -1 00897 index_c = g4x6_list ( iconst ) % c + first_atom -1 00898 index_d = g4x6_list ( iconst ) % d + first_atom -1 00899 v12 ( : ) = vel ( :, index_a ) - vel ( :, index_b ) 00900 v13 ( : ) = vel ( :, index_a ) - vel ( :, index_c ) 00901 v14 ( : ) = vel ( :, index_a ) - vel ( :, index_d ) 00902 v23 ( : ) = vel ( :, index_b ) - vel ( :, index_c ) 00903 v24 ( : ) = vel ( :, index_b ) - vel ( :, index_d ) 00904 v34 ( : ) = vel ( :, index_c ) - vel ( :, index_d ) 00905 00906 r12 ( : ) = particle_set ( index_a ) % r ( : ) - particle_set ( index_b ) % r ( : ) 00907 r13 ( : ) = particle_set ( index_a ) % r ( : ) - particle_set ( index_c ) % r ( : ) 00908 r14 ( : ) = particle_set ( index_a ) % r ( : ) - particle_set ( index_d ) % r ( : ) 00909 r23 ( : ) = particle_set ( index_b ) % r ( : ) - particle_set ( index_c ) % r ( : ) 00910 r24 ( : ) = particle_set ( index_b ) % r ( : ) - particle_set ( index_d ) % r ( : ) 00911 r34 ( : ) = particle_set ( index_c ) % r ( : ) - particle_set ( index_d ) % r ( : ) 00912 atomic_kind => particle_set ( index_a ) % atomic_kind 00913 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 00914 imass1 = 1.0_dp/mass 00915 atomic_kind => particle_set ( index_b ) % atomic_kind 00916 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 00917 imass2 = 1.0_dp/mass 00918 atomic_kind => particle_set ( index_c ) % atomic_kind 00919 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 00920 imass3 = 1.0_dp/mass 00921 atomic_kind => particle_set ( index_d ) % atomic_kind 00922 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 00923 imass4 = 1.0_dp/mass 00924 lg4x6 ( iconst ) % fa = -2.0_dp* r12 00925 lg4x6 ( iconst ) % fb = -2.0_dp* r13 00926 lg4x6 ( iconst ) % fc = -2.0_dp* r14 00927 lg4x6 ( iconst ) % fd = -2.0_dp* r23 00928 lg4x6 ( iconst ) % fe = -2.0_dp* r24 00929 lg4x6 ( iconst ) % ff = -2.0_dp* r34 00930 ! Check for fixed atom constraints 00931 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4,& 00932 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 00933 ! construct matrix 00934 amat(1,1)=(imass1+imass2)*DOTPROD_3D(r12,lg4x6 ( iconst ) % fa ) 00935 amat(1,2)=imass1*DOTPROD_3D(r12, lg4x6 ( iconst) % fb ) 00936 amat(1,3)=imass1*DOTPROD_3D(r12, lg4x6 ( iconst ) % fc ) 00937 amat(1,4)=-imass2*DOTPROD_3D(r12, lg4x6 ( iconst ) % fd ) 00938 amat(1,5)=-imass2*DOTPROD_3D(r12, lg4x6 ( iconst ) % fe ) 00939 amat(1,6)=0.0_dp 00940 00941 amat(2,1)=imass1*DOTPROD_3D(r13,lg4x6 ( iconst ) % fa ) 00942 amat(2,2)=(imass1+imass3)*DOTPROD_3D(r13, lg4x6 ( iconst) % fb ) 00943 amat(2,3)=imass1*DOTPROD_3D(r13, lg4x6 ( iconst ) % fc ) 00944 amat(2,4)=imass3*DOTPROD_3D(r13, lg4x6 ( iconst ) % fd ) 00945 amat(2,5)=0.0_dp 00946 amat(2,6)=-imass3*DOTPROD_3D(r13, lg4x6 ( iconst ) % ff ) 00947 00948 amat(3,1)=imass1*DOTPROD_3D(r14,lg4x6 ( iconst ) % fa ) 00949 amat(3,2)=imass1*DOTPROD_3D(r14, lg4x6 ( iconst) % fb ) 00950 amat(3,3)=(imass1+imass4)*DOTPROD_3D(r14, lg4x6 ( iconst ) % fc ) 00951 amat(3,4)=0.0_dp 00952 amat(3,5)=imass4*DOTPROD_3D(r14, lg4x6 ( iconst ) % fe ) 00953 amat(3,6)=imass4*DOTPROD_3D(r14, lg4x6 ( iconst ) % ff ) 00954 00955 amat(4,1)=-imass2*DOTPROD_3D(r23,lg4x6 ( iconst ) % fa ) 00956 amat(4,2)=imass3*DOTPROD_3D(r23, lg4x6 ( iconst) % fb ) 00957 amat(4,3)=0.0_dp 00958 amat(4,4)=(imass3+imass2)*DOTPROD_3D(r23, lg4x6 ( iconst ) % fd ) 00959 amat(4,5)=imass2*DOTPROD_3D(r23, lg4x6 ( iconst ) % fe ) 00960 amat(4,6)=-imass3*DOTPROD_3D(r23, lg4x6 ( iconst ) % ff ) 00961 00962 amat(5,1)=-imass2*DOTPROD_3D(r24,lg4x6 ( iconst ) % fa ) 00963 amat(5,2)=0.0_dp 00964 amat(5,3)=imass4*DOTPROD_3D(r24, lg4x6 ( iconst) % fc ) 00965 amat(5,4)=imass2*DOTPROD_3D(r24, lg4x6 ( iconst ) % fd ) 00966 amat(5,5)=(imass4+imass2)*DOTPROD_3D(r24, lg4x6 ( iconst ) % fe ) 00967 amat(5,6)=imass4*DOTPROD_3D(r24, lg4x6 ( iconst ) % ff ) 00968 00969 amat(6,1)=0.0_dp 00970 amat(6,2)=-imass3*DOTPROD_3D(r34,lg4x6 ( iconst ) % fb ) 00971 amat(6,3)=imass4*DOTPROD_3D(r34, lg4x6 ( iconst) % fc ) 00972 amat(6,4)=-imass3*DOTPROD_3D(r34, lg4x6 ( iconst ) % fd ) 00973 amat(6,5)=imass4*DOTPROD_3D(r34, lg4x6 ( iconst ) % fe ) 00974 amat(6,6)=(imass3+imass4)*DOTPROD_3D(r34, lg4x6 ( iconst ) % ff ) 00975 00976 ! construct solution vector 00977 bvec(1,1)=DOTPROD_3D(r12,v12) 00978 bvec(2,1)=DOTPROD_3D(r13,v13) 00979 bvec(3,1)=DOTPROD_3D(r14,v14) 00980 bvec(4,1)=DOTPROD_3D(r23,v23) 00981 bvec(5,1)=DOTPROD_3D(r24,v24) 00982 bvec(6,1)=DOTPROD_3D(r34,v34) 00983 bvec=-bvec*2.0_dp*idt 00984 00985 ! get lambda 00986 CALL solve_system ( amat, 6, bvec ) 00987 lg4x6 ( iconst ) % lambda ( : )=bvec ( :, 1 ) 00988 00989 fc1=lg4x6 ( iconst ) % lambda ( 1 ) * lg4x6 ( iconst ) % fa + & 00990 lg4x6 ( iconst ) % lambda ( 2 ) * lg4x6 ( iconst ) % fb + & 00991 lg4x6 ( iconst ) % lambda ( 3 ) * lg4x6 ( iconst ) % fc 00992 fc2=-lg4x6 ( iconst ) % lambda ( 1 ) * lg4x6 ( iconst ) % fa + & 00993 lg4x6 ( iconst ) % lambda ( 4 ) * lg4x6 ( iconst ) % fd + & 00994 lg4x6 ( iconst ) % lambda ( 5 ) * lg4x6 ( iconst ) % fe 00995 fc3=-lg4x6 ( iconst ) % lambda ( 2 ) * lg4x6 ( iconst ) % fb - & 00996 lg4x6 ( iconst ) % lambda ( 4 ) * lg4x6 ( iconst ) % fd + & 00997 lg4x6 ( iconst ) % lambda ( 6 ) * lg4x6 ( iconst ) % ff 00998 fc4=-lg4x6 ( iconst ) % lambda ( 3 ) * lg4x6 ( iconst ) % fc - & 00999 lg4x6 ( iconst ) % lambda ( 5 ) * lg4x6 ( iconst ) % fe - & 01000 lg4x6 ( iconst ) % lambda ( 6 ) * lg4x6 ( iconst ) % ff 01001 vel ( :, index_a ) = vel ( :, index_a ) + imass1 * dtby2 * fc1 ( : ) 01002 vel ( :, index_b ) = vel ( :, index_b ) + imass2 * dtby2 * fc2 ( : ) 01003 vel ( :, index_c ) = vel ( :, index_c ) + imass3 * dtby2 * fc3 ( : ) 01004 vel ( :, index_d ) = vel ( :, index_d ) + imass4 * dtby2 * fc4 ( : ) 01005 END DO 01006 END SUBROUTINE rattle_4x6_low 01007 01008 ! ***************************************************************************** 01012 SUBROUTINE rattle_roll_4x6_low( fixd_list, g4x6_list, lg4x6, first_atom, & 01013 particle_set, vel, r_rattle, dt, veps, error ) 01014 01015 TYPE(fixd_constraint_type), 01016 DIMENSION(:), POINTER :: fixd_list 01017 TYPE(g4x6_constraint_type), POINTER :: g4x6_list( : ) 01018 TYPE(local_g4x6_constraint_type), 01019 POINTER :: lg4x6( : ) 01020 INTEGER, INTENT(IN) :: first_atom 01021 TYPE(particle_type), POINTER :: particle_set( : ) 01022 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 01023 REAL(KIND=dp), DIMENSION(:, :), 01024 INTENT(IN) :: r_rattle 01025 REAL(kind=dp), INTENT(in) :: dt 01026 REAL(KIND=dp), DIMENSION(:, :), 01027 INTENT(IN) :: veps 01028 TYPE(cp_error_type), INTENT(inout) :: error 01029 01030 INTEGER :: iconst, index_a, index_b, 01031 index_c, index_d 01032 REAL(KIND=dp) :: dtby2, idt, imass1, imass2, 01033 imass3, imass4, mass 01034 REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, f_roll4, 01035 f_roll5, f_roll6, fc1, fc2, fc3, fc4, r12, r13, r14, r23, r24, r34, 01036 v12, v13, v14, v23, v24, v34, vec 01037 REAL(KIND=dp), DIMENSION(6) :: lambda 01038 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 01039 REAL(KIND=dp), DIMENSION(6, 6) :: amat 01040 TYPE(atomic_kind_type), POINTER :: atomic_kind 01041 01042 idt = 1.0_dp/dt 01043 dtby2 = dt*.5_dp 01044 DO iconst = 1, SIZE ( g4x6_list ) 01045 IF (g4x6_list(iconst)%restraint%active) CYCLE 01046 index_a = g4x6_list ( iconst ) % a + first_atom -1 01047 index_b = g4x6_list ( iconst ) % b + first_atom -1 01048 index_c = g4x6_list ( iconst ) % c + first_atom -1 01049 index_d = g4x6_list ( iconst ) % d + first_atom -1 01050 v12 ( : ) = vel ( :, index_a ) - vel ( :, index_b ) 01051 v13 ( : ) = vel ( :, index_a ) - vel ( :, index_c ) 01052 v14 ( : ) = vel ( :, index_a ) - vel ( :, index_d ) 01053 v23 ( : ) = vel ( :, index_b ) - vel ( :, index_c ) 01054 v24 ( : ) = vel ( :, index_b ) - vel ( :, index_d ) 01055 v34 ( : ) = vel ( :, index_c ) - vel ( :, index_d ) 01056 01057 r12 ( : ) = particle_set ( index_a ) % r ( : ) - particle_set ( index_b ) % r ( : ) 01058 r13 ( : ) = particle_set ( index_a ) % r ( : ) - particle_set ( index_c ) % r ( : ) 01059 r14 ( : ) = particle_set ( index_a ) % r ( : ) - particle_set ( index_d ) % r ( : ) 01060 r23 ( : ) = particle_set ( index_b ) % r ( : ) - particle_set ( index_c ) % r ( : ) 01061 r24 ( : ) = particle_set ( index_b ) % r ( : ) - particle_set ( index_d ) % r ( : ) 01062 r34 ( : ) = particle_set ( index_c ) % r ( : ) - particle_set ( index_d ) % r ( : ) 01063 atomic_kind => particle_set ( index_a ) % atomic_kind 01064 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 01065 imass1 = 1.0_dp/mass 01066 atomic_kind => particle_set ( index_b ) % atomic_kind 01067 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 01068 imass2 = 1.0_dp/mass 01069 atomic_kind => particle_set ( index_c ) % atomic_kind 01070 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 01071 imass3 = 1.0_dp/mass 01072 atomic_kind => particle_set ( index_d ) % atomic_kind 01073 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass) 01074 imass4 = 1.0_dp/mass 01075 lg4x6 ( iconst ) % fa = -2.0_dp* r12 01076 lg4x6 ( iconst ) % fb = -2.0_dp* r13 01077 lg4x6 ( iconst ) % fc = -2.0_dp* r14 01078 lg4x6 ( iconst ) % fd = -2.0_dp* r23 01079 lg4x6 ( iconst ) % fe = -2.0_dp* r24 01080 lg4x6 ( iconst ) % ff = -2.0_dp* r34 01081 ! Check for fixed atom constraints 01082 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4,& 01083 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 01084 ! roll the fc 01085 CALL MATVEC_3x3 (f_roll1, r_rattle, lg4x6 ( iconst ) % fa ) 01086 CALL MATVEC_3x3 (f_roll2, r_rattle, lg4x6 ( iconst ) % fb ) 01087 CALL MATVEC_3x3 (f_roll3, r_rattle, lg4x6 ( iconst ) % fc ) 01088 CALL MATVEC_3x3 (f_roll4, r_rattle, lg4x6 ( iconst ) % fd ) 01089 CALL MATVEC_3x3 (f_roll5, r_rattle, lg4x6 ( iconst ) % fe ) 01090 CALL MATVEC_3x3 (f_roll6, r_rattle, lg4x6 ( iconst ) % ff ) 01091 ! construct matrix 01092 amat(1,1)=(imass1+imass2)*DOTPROD_3D(r12,f_roll1 ) 01093 amat(1,2)=imass1*DOTPROD_3D(r12, f_roll2 ) 01094 amat(1,3)=imass1*DOTPROD_3D(r12, f_roll3 ) 01095 amat(1,4)=-imass2*DOTPROD_3D(r12, f_roll4 ) 01096 amat(1,5)=-imass2*DOTPROD_3D(r12, f_roll5 ) 01097 amat(1,6)=0.0_dp 01098 01099 amat(2,1)=imass1*DOTPROD_3D(r13,f_roll1 ) 01100 amat(2,2)=(imass1+imass3)*DOTPROD_3D(r13, f_roll2) 01101 amat(2,3)=imass1*DOTPROD_3D(r13, f_roll3 ) 01102 amat(2,4)=imass3*DOTPROD_3D(r13, f_roll4 ) 01103 amat(2,5)=0.0_dp 01104 amat(2,6)=-imass3*DOTPROD_3D(r13, f_roll6 ) 01105 01106 amat(3,1)=imass1*DOTPROD_3D(r14,f_roll1 ) 01107 amat(3,2)=imass1*DOTPROD_3D(r14, f_roll2 ) 01108 amat(3,3)=(imass1+imass4)*DOTPROD_3D(r14, f_roll3 ) 01109 amat(3,4)=0.0_dp 01110 amat(3,5)=imass4*DOTPROD_3D(r14, f_roll5 ) 01111 amat(3,6)=imass4*DOTPROD_3D(r14, f_roll6 ) 01112 01113 amat(4,1)=-imass2*DOTPROD_3D(r23,f_roll1 ) 01114 amat(4,2)=imass3*DOTPROD_3D(r23, f_roll2 ) 01115 amat(4,3)=0.0_dp 01116 amat(4,4)=(imass3+imass2)*DOTPROD_3D(r23, f_roll4 ) 01117 amat(4,5)=imass2*DOTPROD_3D(r23, f_roll5 ) 01118 amat(4,6)=-imass3*DOTPROD_3D(r23, f_roll6 ) 01119 01120 amat(5,1)=-imass2*DOTPROD_3D(r24,f_roll1 ) 01121 amat(5,2)=0.0_dp 01122 amat(5,3)=imass4*DOTPROD_3D(r24, f_roll3) 01123 amat(5,4)=imass2*DOTPROD_3D(r24, f_roll4 ) 01124 amat(5,5)=(imass4+imass2)*DOTPROD_3D(r24, f_roll5 ) 01125 amat(5,6)=imass4*DOTPROD_3D(r24, f_roll6 ) 01126 01127 amat(6,1)=0.0_dp 01128 amat(6,2)=-imass3*DOTPROD_3D(r34,f_roll2 ) 01129 amat(6,3)=imass4*DOTPROD_3D(r34, f_roll3 ) 01130 amat(6,4)=-imass3*DOTPROD_3D(r34, f_roll4 ) 01131 amat(6,5)=imass4*DOTPROD_3D(r34, f_roll5 ) 01132 amat(6,6)=(imass3+imass4)*DOTPROD_3D(r34, f_roll6 ) 01133 01134 ! construct solution vector 01135 CALL matvec_3x3 (vec, veps, r12 ) 01136 bvec ( 1, 1 ) = DOTPROD_3D ( r12, v12 + vec ) 01137 CALL matvec_3x3 (vec, veps, r13 ) 01138 bvec ( 2, 1 ) = DOTPROD_3D ( r13, v13 + vec ) 01139 CALL matvec_3x3 (vec, veps, r14 ) 01140 bvec ( 3, 1 ) = DOTPROD_3D ( r14, v14 + vec ) 01141 CALL matvec_3x3 (vec, veps, r23 ) 01142 bvec ( 4, 1 ) = DOTPROD_3D ( r23, v23 + vec ) 01143 CALL matvec_3x3 (vec, veps, r24 ) 01144 bvec ( 5, 1 ) = DOTPROD_3D ( r24, v24 + vec ) 01145 CALL matvec_3x3 (vec, veps, r34 ) 01146 bvec ( 6, 1 ) = DOTPROD_3D ( r34, v34 + vec ) 01147 bvec=-bvec*2.0_dp*idt 01148 01149 ! get lambda 01150 CALL solve_system ( amat, 6, bvec ) 01151 lambda ( : ) = bvec ( :, 1 ) 01152 lg4x6 ( iconst ) % lambda ( : )= lambda 01153 01154 fc1=lambda ( 1 ) * f_roll1 + & 01155 lambda ( 2 ) * f_roll2 + & 01156 lambda ( 3 ) * f_roll3 01157 fc2=-lambda ( 1 ) * f_roll1 + & 01158 lambda ( 4 ) * f_roll4 + & 01159 lambda ( 5 ) * f_roll5 01160 fc3=-lambda ( 2 ) * f_roll2 - & 01161 lambda ( 4 ) * f_roll4 + & 01162 lambda ( 6 ) * f_roll6 01163 fc4=-lambda ( 3 ) * f_roll3 - & 01164 lambda ( 5 ) * f_roll5 - & 01165 lambda ( 6 ) * f_roll6 01166 vel ( :, index_a ) = vel ( :, index_a ) + imass1 * dtby2 * fc1 ( : ) 01167 vel ( :, index_b ) = vel ( :, index_b ) + imass2 * dtby2 * fc2 ( : ) 01168 vel ( :, index_c ) = vel ( :, index_c ) + imass3 * dtby2 * fc3 ( : ) 01169 vel ( :, index_d ) = vel ( :, index_d ) + imass4 * dtby2 * fc4 ( : ) 01170 END DO 01171 END SUBROUTINE rattle_roll_4x6_low 01172 01173 END MODULE constraint_4x6
1.7.3