hdiff output

r33130/chirality.F90 2017-08-07 15:30:55.513589445 +0100 r33129/chirality.F90 2017-08-07 15:30:57.057609985 +0100
195:    integer                                      :: file_length195:    integer                                      :: file_length
196: 196: 
197: #ifdef _SVN_ROOT_197: #ifdef _SVN_ROOT_
198:    call system('python ' // _SVN_ROOT_ // '/SCRIPTS/AMBER/chirality/chirality.py' // ' coords.prmtop')198:    call system('python ' // _SVN_ROOT_ // '/SCRIPTS/AMBER/chirality/chirality.py' // ' coords.prmtop')
199: #else199: #else
200:    call system('python ' // chirality_script // ' coords.prmtop')200:    call system('python ' // chirality_script // ' coords.prmtop')
201: #endif201: #endif
202: 202: 
203: ! Work out the number of chiral centres by reading the .chirality_list file 203: ! Work out the number of chiral centres by reading the .chirality_list file 
204:    num_chiral_centres = file_length('.chirality_list')204:    num_chiral_centres = file_length('.chirality_list')
205:    if (.not. allocated(sr_atoms)) allocate(sr_atoms(num_chiral_centres, 5))205:    if (allocated(sr_atoms)) deallocate(sr_atoms)
 206:    allocate(sr_atoms(num_chiral_centres, 5))
206: 207: 
207: ! Now read the chiral centres into sr_atoms208: ! Now read the chiral centres into sr_atoms
208:    call file_open('.chirality_list', file_unit, .false.)209:    call file_open('.chirality_list', file_unit, .false.)
209:    do i = 1, num_chiral_centres210:    do i = 1, num_chiral_centres
210:       read(file_unit, '(5i8)') sr_atoms(i, :)211:       read(file_unit, '(5i8)') sr_atoms(i, :)
211:    end do212:    end do
212:    close(file_unit)213:    close(file_unit)
213: 214: 
214: ! Print a test copy215: ! Print a test copy
215: !   call file_open('chirality_list_copy', file_unit, .false.)216: !   call file_open('chirality_list_copy', file_unit, .false.)
216: !   do i = 1, num_chiral_centres217: !   do i = 1, num_chiral_centres
217: !      write(file_unit, '(5i10)') sr_atoms(i, :)218: !      write(file_unit, '(5i10)') sr_atoms(i, :)
218: !   end do219: !   end do
219: !   close(file_unit)220: !   close(file_unit)
220: 221: 
221: ! Now calculate the chirality of the centres and save it in sr_states_initial222: ! Now calculate the chirality of the centres and save it in sr_states_initial
222:    if (.not. allocated(sr_states_initial)) allocate(sr_states_initial(num_chiral_centres))223:    if (allocated(sr_states_initial)) deallocate(sr_states_initial)
 224:    allocate(sr_states_initial(num_chiral_centres))
223:    do i = 1, num_chiral_centres225:    do i = 1, num_chiral_centres
224:       atom_number = sr_atoms(i, 1)226:       atom_number = sr_atoms(i, 1)
225:       centre_coords(1) = coords(3 * atom_number - 2)227:       centre_coords(1) = coords(3 * atom_number - 2)
226:       centre_coords(2) = coords(3 * atom_number - 1)228:       centre_coords(2) = coords(3 * atom_number - 1)
227:       centre_coords(3) = coords(3 * atom_number    )229:       centre_coords(3) = coords(3 * atom_number    )
228:       do j = 1, 4230:       do j = 1, 4
229:          atom_number = sr_atoms(i, j + 1) 231:          atom_number = sr_atoms(i, j + 1) 
230:          neighbour_coords(3 * j - 2) = coords(3 * atom_number - 2)232:          neighbour_coords(3 * j - 2) = coords(3 * atom_number - 2)
231:          neighbour_coords(3 * j - 1) = coords(3 * atom_number - 1)233:          neighbour_coords(3 * j - 1) = coords(3 * atom_number - 1)
232:          neighbour_coords(3 * j    ) = coords(3 * atom_number    )234:          neighbour_coords(3 * j    ) = coords(3 * atom_number    )
269:    integer                                      :: file_length271:    integer                                      :: file_length
270: 272: 
271: #ifdef _SVN_ROOT_273: #ifdef _SVN_ROOT_
272:    call system('python ' // _SVN_ROOT_ // '/SCRIPTS/AMBER/chirality/cistrans.py' // ' coords.prmtop')274:    call system('python ' // _SVN_ROOT_ // '/SCRIPTS/AMBER/chirality/cistrans.py' // ' coords.prmtop')
273: #else275: #else
274:    call system('python ' // cis_trans_script // ' coords.prmtop')276:    call system('python ' // cis_trans_script // ' coords.prmtop')
275: #endif277: #endif
276: 278: 
277: ! Work out the number of peptide bonds by reading the .cis_trans_list file 279: ! Work out the number of peptide bonds by reading the .cis_trans_list file 
278:    num_peptide_bonds = file_length('.cis_trans_list')280:    num_peptide_bonds = file_length('.cis_trans_list')
279:    if (.not. allocated(cis_trans_atoms)) allocate(cis_trans_atoms(num_peptide_bonds, 4))281:    if (allocated(cis_trans_atoms)) deallocate(cis_trans_atoms)
 282:    allocate(cis_trans_atoms(num_peptide_bonds, 4))
280: 283: 
281: ! Now read the chiral centres into cis_trans_atoms284: ! Now read the chiral centres into cis_trans_atoms
282:    call file_open('.cis_trans_list', file_unit, .false.)285:    call file_open('.cis_trans_list', file_unit, .false.)
283:    do i = 1, num_peptide_bonds286:    do i = 1, num_peptide_bonds
284:       read(file_unit, '(4i8)') cis_trans_atoms(i, :)287:       read(file_unit, '(4i8)') cis_trans_atoms(i, :)
285:    end do288:    end do
286:    close(file_unit)289:    close(file_unit)
287: 290: 
288: ! Print a test copy291: ! Print a test copy
289: !   call file_open('cis_trans_list_copy', file_unit, .false.)292: !   call file_open('cis_trans_list_copy', file_unit, .false.)
290: !   do i = 1, num_peptide_bonds293: !   do i = 1, num_peptide_bonds
291: !      write(file_unit, '(4i8)') cis_trans_atoms(i, :)294: !      write(file_unit, '(4i8)') cis_trans_atoms(i, :)
292: !   end do295: !   end do
293: !   close(file_unit)296: !   close(file_unit)
294: 297: 
295: ! Now calculate the isomerism of the peptide bonds and save it in cis_trans_states_initial298: ! Now calculate the isomerism of the peptide bonds and save it in cis_trans_states_initial
296:    if (.not. allocated(cis_trans_states_initial)) allocate(cis_trans_states_initial(num_peptide_bonds))299:    if (allocated(cis_trans_states_initial)) deallocate(cis_trans_states_initial)
 300:    allocate(cis_trans_states_initial(num_peptide_bonds))
297:    do i = 1, num_peptide_bonds301:    do i = 1, num_peptide_bonds
298:       do j = 1, 4302:       do j = 1, 4
299:          atom_number = cis_trans_atoms(i, j) 303:          atom_number = cis_trans_atoms(i, j) 
300:          peptide_coords(3 * j - 2) = coords(3 * atom_number - 2)304:          peptide_coords(3 * j - 2) = coords(3 * atom_number - 2)
301:          peptide_coords(3 * j - 1) = coords(3 * atom_number - 1)305:          peptide_coords(3 * j - 1) = coords(3 * atom_number - 1)
302:          peptide_coords(3 * j    ) = coords(3 * atom_number    )306:          peptide_coords(3 * j    ) = coords(3 * atom_number    )
303:       end do307:       end do
304:       cis_trans_states_initial(i) = cis_trans(peptide_coords)308:       cis_trans_states_initial(i) = cis_trans(peptide_coords)
305:    end do 309:    end do 
306: 310: 
544: 548: 
545:    print *, "Angle:" , 180.0 * dihedral(plane_coords_1) / pi549:    print *, "Angle:" , 180.0 * dihedral(plane_coords_1) / pi
546:    print *, "Angle:" , 180.0 * dihedral(plane_coords_2) / pi550:    print *, "Angle:" , 180.0 * dihedral(plane_coords_2) / pi
547:    print *, "Angle:", 180.0 * dihedral(right_angle_coords_1) / pi551:    print *, "Angle:", 180.0 * dihedral(right_angle_coords_1) / pi
548:    print *, "Angle:", 180.0 * dihedral(right_angle_coords_2) / pi552:    print *, "Angle:", 180.0 * dihedral(right_angle_coords_2) / pi
549: 553: 
550:    print *, "Right angle cis/trans: ", cis_trans(right_angle_coords_1)554:    print *, "Right angle cis/trans: ", cis_trans(right_angle_coords_1)
551: 555: 
552: end subroutine test_chirality556: end subroutine test_chirality
553: 557: 
 558: subroutine dealloc_states_mutation()
 559:    if (allocated(sr_states)) deallocate(sr_states)
 560: end subroutine dealloc_states_mutation
 561: 
554: end module chirality562: end module chirality


r33130/gay-berne.f90 2017-08-07 15:30:55.741592478 +0100 r33129/gay-berne.f90 2017-08-07 15:30:57.285613019 +0100
5292:       MG = (1.D0 - LAMDA) * AEINV + LAMDA * BEINV5292:       MG = (1.D0 - LAMDA) * AEINV + LAMDA * BEINV
5293: 5293: 
5294:       CALL MTRXIN(MG, MGINV)5294:       CALL MTRXIN(MG, MGINV)
5295: 5295: 
5296:       MGINVR  =  MATMUL(MGINV, RIJ) 5296:       MGINVR  =  MATMUL(MGINV, RIJ) 
5297: 5297: 
5298:       SLMD =  - LAMDA * (1.D0 - LAMDA) * DOT_PRODUCT(RIJ,MGINVR)5298:       SLMD =  - LAMDA * (1.D0 - LAMDA) * DOT_PRODUCT(RIJ,MGINVR)
5299: 5299: 
5300:       RETURN5300:       RETURN
5301:       END SUBROUTINE OBJCTF5301:       END SUBROUTINE OBJCTF
 5302: 
 5303: SUBROUTINE GENERATECGDIMER
 5304: ! this will be useful if we have a reference dimer and would like to get the
 5305: ! coarse-grained coordinates for it (virus capsid protein dimers, for example),
 5306: ! by starting from the angle-axis information about the rigid monomer containing
 5307: ! an arbitrary number of ellipsoids. For this, we need as input the coordinate centres for both units,
 5308: ! and the rotation matrix describing the rotation of the first unit to the second one. Output 
 5309: ! will be the angle-axis coordinates of the second unit.
 5310: ! Individual orientations of ellipsoids should go in the pysites.xyz file.
 5311: use rotations
 5312: implicit none
 5313: integer i,j,k, realnatoms, nlines
 5314: double precision :: y(12),atomblockcentre(2,3), rotationmatrix(3,3)
 5315: 
 5316: open(255,file='coords',status='old')
 5317: 
 5318: call determinelines(255,nlines)
 5319: if(nlines>4) write(*,*) 'ERROR - we should be generating a rigid body dimer, so we need 4 coordinates in the coords file'
 5320: 
 5321: do i=1,nlines
 5322:    read(255,*) y(3*i-2), y(3*i-1), y(3*i)
 5323: end do
 5324: 
 5325: close(255)
 5326: 
 5327: open(256,file='rotationmatrix',status='old')
 5328:  do j=1,3
 5329:    read(256,*) rotationmatrix(1:3,j)
 5330:  end do
 5331:  write(*,*) 'read rotation matrix:'
 5332:  write(*,*) rotationmatrix(:,:)
 5333:  realnatoms=nlines/2
 5334:     y(4:6)=MATMUL(y(1:3),rotationmatrix)+y(4:6) 
 5335:     y(3*2+3*REALNATOMS-2:3*2+3*REALNATOMS)=rot_rotate_aa(y(3*2+3*REALNATOMS-2:3*2+3*REALNATOMS),rot_mx2aa(rotationmatrix(:,:)))
 5336: write(*,*) 'finished generating dimer coordinates'
 5337: do i=1,nlines
 5338:    write(*,*) y(3*i-2), y(3*i-1), y(3*i)
 5339: end do
 5340: 
 5341: contains
 5342: ! Returns the inverse of a matrix calculated by finding the LU
 5343: ! decomposition.  Depends on LAPACK.
 5344: function inv(A) result(Ainv)
 5345:   double precision, dimension(:,:), intent(in) :: A
 5346:   double precision, dimension(size(A,1),size(A,2)) :: Ainv
 5347: 
 5348:   double precision, dimension(size(A,1)) :: work  ! work array for LAPACK
 5349:   integer, dimension(size(A,1)) :: ipiv   ! pivot indices
 5350:   integer :: n, info
 5351: 
 5352:   ! External procedures defined in LAPACK
 5353:   external DGETRF
 5354:   external DGETRI
 5355: 
 5356:   ! Store A in Ainv to prevent it from being overwritten by LAPACK
 5357:   Ainv = A
 5358:   n = size(A,1)
 5359: 
 5360:   ! DGETRF computes an LU factorization of a general M-by-N matrix A
 5361:   ! using partial pivoting with row interchanges.
 5362:   call DGETRF(n, n, Ainv, n, ipiv, info)
 5363: 
 5364:   if (info /= 0) then
 5365:      stop 'Matrix is numerically singular!'
 5366:   end if
 5367: 
 5368:   ! DGETRI computes the inverse of a matrix using the LU factorization
 5369:   ! computed by DGETRF.
 5370:   call DGETRI(n, Ainv, n, ipiv, work, n, info)
 5371: 
 5372:   if (info /= 0) then
 5373:      stop 'Matrix inversion failed!'
 5374:   end if
 5375: end function inv
 5376: 
 5377: 
 5378: END SUBROUTINE GENERATECGDIMER
 5379: 
 5380: 
 5381: SUBROUTINE DETERMINELINES(nunit,nlines)
 5382: implicit none
 5383: integer nunit, nlines, iostatus
 5384: character(len=10) check
 5385: 
 5386: REWIND(nunit)
 5387: 
 5388: nlines=0
 5389: do
 5390:   IF(iostatus<0) EXIT
 5391:   nlines = nlines + 1
 5392:   READ(nunit,*,iostat=iostatus) check
 5393: !  write(*,*) check,nunit
 5394: end do
 5395:   nlines = nlines - 1
 5396:   REWIND(nunit)
 5397: RETURN
 5398: 
 5399: 
 5400: END SUBROUTINE DETERMINELINES
 5401: 
 5402: 


r33130/io1.f 2017-08-07 15:30:55.961595405 +0100 r33129/io1.f 2017-08-07 15:30:57.617617444 +0100
115:                    READ(COORDS_UNIT,*) COORDS(J1,JP)115:                    READ(COORDS_UNIT,*) COORDS(J1,JP)
116:                ENDDO116:                ENDDO
117:             ENDDO117:             ENDDO
118:          ELSEIF (MLQT) THEN ! for this ML it is one variable per line118:          ELSEIF (MLQT) THEN ! for this ML it is one variable per line
119:             REWIND(COORDS_UNIT)119:             REWIND(COORDS_UNIT)
120:             DO JP=1,NPAR120:             DO JP=1,NPAR
121:                DO J1=1,NATOMS121:                DO J1=1,NATOMS
122:                    READ(COORDS_UNIT,*) COORDS(J1,JP)122:                    READ(COORDS_UNIT,*) COORDS(J1,JP)
123:                ENDDO123:                ENDDO
124:             ENDDO124:             ENDDO
 125:          ELSEIF (ORBITALS) THEN ! for this orbital landscape it is one variable per line
 126:             REWIND(COORDS_UNIT)
 127:             DO JP=1,NPAR
 128:                DO J1=1,NATOMS
 129:                    READ(COORDS_UNIT,*) COORDS(J1,JP)
 130:                ENDDO
 131:             ENDDO
125:          ELSE132:          ELSE
126:             REWIND(COORDS_UNIT)133:             REWIND(COORDS_UNIT)
127:             DO JP=1,NPAR134:             DO JP=1,NPAR
128:                DO J1=1,NATOMS135:                DO J1=1,NATOMS
129:                   J2=3*(J1-1)136:                   J2=3*(J1-1)
130:                    READ(COORDS_UNIT,*) COORDS(J2+1,JP), COORDS(J2+2,JP), COORDS(J2+3,JP)137:                    READ(COORDS_UNIT,*) COORDS(J2+1,JP), COORDS(J2+2,JP), COORDS(J2+3,JP)
131:                ENDDO138:                ENDDO
132:             ENDDO139:             ENDDO
133:          ENDIF140:          ENDIF
134:          CLOSE(COORDS_UNIT)141:          CLOSE(COORDS_UNIT)
195: 202: 
196:       IF (.NOT.SEEDT.AND..NOT.AMHT.AND..NOT.SUPPRESST) THEN203:       IF (.NOT.SEEDT.AND..NOT.AMHT.AND..NOT.SUPPRESST) THEN
197:          WRITE(MYUNIT,20) 204:          WRITE(MYUNIT,20) 
198: 20       FORMAT('Initial coordinates:')205: 20       FORMAT('Initial coordinates:')
199:          IF (MPIT) THEN206:          IF (MPIT) THEN
200:             WRITE(MYUNIT,30) (COORDS(J1,MYNODE+1),J1=1,3*NATOMS)207:             WRITE(MYUNIT,30) (COORDS(J1,MYNODE+1),J1=1,3*NATOMS)
201:          ELSEIF (MLP3T.OR.MLPVB3T) THEN 208:          ELSEIF (MLP3T.OR.MLPVB3T) THEN 
202:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)209:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)
203:          ELSEIF (MLQT) THEN 210:          ELSEIF (MLQT) THEN 
204:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)211:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)
 212:          ELSEIF (ORBITALS) THEN 
 213:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)
205:          ELSE 214:          ELSE 
206:            DO JP=1,NPAR215:            DO JP=1,NPAR
207:                WRITE(MYUNIT,30) (COORDS(J1,JP),J1=1,3*NATOMS)216:                WRITE(MYUNIT,30) (COORDS(J1,JP),J1=1,3*NATOMS)
208: 30             FORMAT(3F20.10)217: 30             FORMAT(3F20.10)
209:             ENDDO218:             ENDDO
210:          ENDIF219:          ENDIF
211:       ENDIF220:       ENDIF
212: 221: 
213:       IF (MSORIGT) THEN222:       IF (MSORIGT) THEN
214:          WRITE(MYUNIT,'(I4,A)') NATOMS,' M and S silicon atoms'223:          WRITE(MYUNIT,'(I4,A)') NATOMS,' M and S silicon atoms'
714:       ELSEIF (MKTRAPT) THEN723:       ELSEIF (MKTRAPT) THEN
715:          WRITE(MYUNIT,'(I4,A)') NATOMS,' MK trapped ions'724:          WRITE(MYUNIT,'(I4,A)') NATOMS,' MK trapped ions'
716:       ELSEIF (DJWRBT) THEN725:       ELSEIF (DJWRBT) THEN
717:          IF (DJWRBID.EQ.1) THEN726:          IF (DJWRBID.EQ.1) THEN
718:             WRITE(MYUNIT,'(3(I4,A))') NATOMS,' sites for ',NRIGIDBODY,' rigid bodies - DJW potential ',DJWRBID727:             WRITE(MYUNIT,'(3(I4,A))') NATOMS,' sites for ',NRIGIDBODY,' rigid bodies - DJW potential ',DJWRBID
719:             WRITE(MYUNIT,'(2(I4,A))') NRIGIDBODY-NHEXAMERS,' pentamers and ',NHEXAMERS,' hexamers'728:             WRITE(MYUNIT,'(2(I4,A))') NRIGIDBODY-NHEXAMERS,' pentamers and ',NHEXAMERS,' hexamers'
720:             WRITE(MYUNIT,'(A,4G20.10)') 'rho, eps, sigma and pentamer radius=',CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT729:             WRITE(MYUNIT,'(A,4G20.10)') 'rho, eps, sigma and pentamer radius=',CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT
721:             WRITE(MYUNIT,'(A,4G20.10)') 'hexamer sigma, radius and hex/pent sigma=',SIGMAHEX,RADHEX,SIGMAPH730:             WRITE(MYUNIT,'(A,4G20.10)') 'hexamer sigma, radius and hex/pent sigma=',SIGMAHEX,RADHEX,SIGMAPH
722:          ENDIF731:          ENDIF
723:       ELSEIF (MLPVB3T) THEN732:       ELSEIF (MLPVB3T) THEN
 733:          IF (MLPVB3NNT) THEN
 734:             DO J1=1,NATOMS
 735:                IF (FROZEN(J1)) COORDS(J1,1)=0.0D0 ! for frozen zero nearest-neighbour setup
 736:             ENDDO
 737:          ENDIF
724:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLPVB3'738:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLPVB3'
725:       ELSEIF (MLP3T) THEN739:       ELSEIF (MLP3T) THEN
726:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLP3'740:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLP3'
727:       ELSEIF (MLQT) THEN741:       ELSEIF (MLQT) THEN
728:          WRITE(MYUNIT,'(I4,A)') NATOMS,' variables for ML quadratic'742:          WRITE(MYUNIT,'(I4,A)') NATOMS,' variables for ML quadratic'
 743:       ELSEIF (ORBITALS) THEN
 744:          WRITE(MYUNIT,'(I4,A)') NATOMS,' rotations for orbital landscape'
729:       ELSEIF  (LJADDT) THEN745:       ELSEIF  (LJADDT) THEN
730:          IF (SORTT) THEN746:          IF (SORTT) THEN
731:             WRITE(MYUNIT,'(A)') 'Turning off SORT option for LJADD'747:             WRITE(MYUNIT,'(A)') 'Turning off SORT option for LJADD'
732:             SORTT=.FALSE.748:             SORTT=.FALSE.
733:          ENDIF749:          ENDIF
734:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ addressable atoms'750:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ addressable atoms'
735:       ELSE751:       ELSE
736:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ atoms'752:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ atoms'
737:       ENDIF753:       ENDIF
738:       IF (PYGPERIODICT.OR.PYBINARYT) CALL INITIALISEPYGPERIODIC754:       IF (PYGPERIODICT.OR.PYBINARYT) CALL INITIALISEPYGPERIODIC
753:          WRITE(MYUNIT, '(A,I6,A)') 'Searching for approximate symmetry elements every ',NSYMINTERVAL,' steps'769:          WRITE(MYUNIT, '(A,I6,A)') 'Searching for approximate symmetry elements every ',NSYMINTERVAL,' steps'
754:          WRITE(MYUNIT, '(A,5F15.5)') 'Distance tolerances: ',SYMTOL1,SYMTOL2,SYMTOL3,SYMTOL4,SYMTOL5770:          WRITE(MYUNIT, '(A,5F15.5)') 'Distance tolerances: ',SYMTOL1,SYMTOL2,SYMTOL3,SYMTOL4,SYMTOL5
755:          WRITE(MYUNIT, '(A,F15.5)') 'Threshold for distinguishing transformation matrices: ',MATDIFF771:          WRITE(MYUNIT, '(A,F15.5)') 'Threshold for distinguishing transformation matrices: ',MATDIFF
756:          WRITE(MYUNIT,'(A,F15.5)') 'Exponential factor in core-weighted centre of mass calculation: ',DISTFAC772:          WRITE(MYUNIT,'(A,F15.5)') 'Exponential factor in core-weighted centre of mass calculation: ',DISTFAC
757:          WRITE(MYUNIT, '(A,I5)') 'Maximum number of quenches for floater/hole permutations=',NSYMQMAX773:          WRITE(MYUNIT, '(A,I5)') 'Maximum number of quenches for floater/hole permutations=',NSYMQMAX
758:          IF (MOVESHELLT) WRITE(MYUNIT,'(A,I8,A,F12.5)') 'Shell moves allowed in blocks of ',SHELLMOVEMAX,' with probability ',774:          IF (MOVESHELLT) WRITE(MYUNIT,'(A,I8,A,F12.5)') 'Shell moves allowed in blocks of ',SHELLMOVEMAX,' with probability ',
759:      &                        SHELLPROB775:      &                        SHELLPROB
760:       ENDIF776:       ENDIF
761:       IF (DEBUG.OR.CHECKMARKOVT) WRITE(MYUNIT,'(A,I6,A)') 'io1> checking the energy of the saved coordinates in the chain'777:       IF (DEBUG.OR.CHECKMARKOVT) WRITE(MYUNIT,'(A,I6,A)') 'io1> checking the energy of the saved coordinates in the chain'
762:       IF (FREEZE) THEN778:       IF (FREEZE) THEN
763:          IF (MLQT.OR.MLPVB3T.OR.MLP3T) THEN779:          IF (MLQT.OR.MLPVB3T.OR.MLP3T.OR.ORBITALS) THEN
764:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' variables will be frozen:'780:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' variables will be frozen:'
765:          ELSE781:          ELSE
766:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' atoms will be frozen:'782:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' atoms will be frozen:'
767:          ENDIF783:          ENDIF
768:          DO J1=1,NATOMS784:          DO J1=1,NATOMS
769:             IF (FROZEN(J1)) WRITE(MYUNIT,'(I6)') J1785:             IF (FROZEN(J1)) WRITE(MYUNIT,'(I6)') J1
770:          ENDDO786:          ENDDO
771:       ENDIF787:       ENDIF
772:       IF (HARMONICF) THEN788:       IF (HARMONICF) THEN
773:          WRITE(MYUNIT,'(A,F12.4)') 'io1> harmonically constrained atoms: strength = ', HARMONICSTR789:          WRITE(MYUNIT,'(A,F12.4)') 'io1> harmonically constrained atoms: strength = ', HARMONICSTR
902: !            WRITE(MYUNIT, '(A,G15.5,A,G15.5)') 'Maximum step size scaled by estimated nearest neighbour distance of ',918: !            WRITE(MYUNIT, '(A,G15.5,A,G15.5)') 'Maximum step size scaled by estimated nearest neighbour distance of ',
903: !    &                    0.677441D0-0.0037582*NATOMS+9.40318D-6*NATOMS**2-6.21931D-9*NATOMS**3,' to give ',STEP(1)919: !    &                    0.677441D0-0.0037582*NATOMS+9.40318D-6*NATOMS**2-6.21931D-9*NATOMS**3,' to give ',STEP(1)
904:          ELSEIF (MULLERBROWNT) THEN 920:          ELSEIF (MULLERBROWNT) THEN 
905:             RADIUS=100.0D0921:             RADIUS=100.0D0
906:          ELSE 922:          ELSE 
907:             RADIUS=RADIUS*2.0D0**(1.0D0/6.0D0)923:             RADIUS=RADIUS*2.0D0**(1.0D0/6.0D0)
908:          ENDIF924:          ENDIF
909:       ENDIF925:       ENDIF
910:       IF ((.NOT.PERIODIC).AND.(.NOT.AMBER).AND.(.NOT.BLNT).AND.(.NOT.MULLERBROWNT).AND.(.NOT.MODEL1T).AND.(.NOT.PERCOLATET) 926:       IF ((.NOT.PERIODIC).AND.(.NOT.AMBER).AND.(.NOT.BLNT).AND.(.NOT.MULLERBROWNT).AND.(.NOT.MODEL1T).AND.(.NOT.PERCOLATET) 
911:      &                    .AND.(.NOT.QCIPOTT).AND.(.NOT.INTCONSTRAINTT).AND.(.NOT.MLP3T).AND.(.NOT.MKTRAPT).AND.(.NOT.MLQT)927:      &                    .AND.(.NOT.QCIPOTT).AND.(.NOT.INTCONSTRAINTT).AND.(.NOT.MLP3T).AND.(.NOT.MKTRAPT).AND.(.NOT.MLQT)
912:      &                    .AND.(.NOT.MLPVB3T)) 928:      &                    .AND.(.NOT.MLPVB3T).AND.(.NOT.ORBITALS)) 
913:      1                    WRITE(MYUNIT,'(A,F20.10)') 'Container radius=',RADIUS929:      1                    WRITE(MYUNIT,'(A,F20.10)') 'Container radius=',RADIUS
914:       RADIUS=RADIUS**2930:       RADIUS=RADIUS**2
915:       IF (PERCOLATET) WRITE(MYUNIT,'(A,F20.10)') 'Checking for percolated structure, cutoff=',PERCCUT931:       IF (PERCOLATET) WRITE(MYUNIT,'(A,F20.10)') 'Checking for percolated structure, cutoff=',PERCCUT
916:       PERCCUT=PERCCUT**2932:       PERCCUT=PERCCUT**2
917:       IF (NPAR.GT.1) THEN933:       IF (NPAR.GT.1) THEN
918:          WRITE(MYUNIT,'(I2,A)') NPAR,' parallel runs'934:          WRITE(MYUNIT,'(I2,A)') NPAR,' parallel runs'
919:          IF (TABOOT) WRITE(MYUNIT,'(A,I4,A)') 'Taboo lists contain the lowest ',NTAB,' minima'935:          IF (TABOOT) WRITE(MYUNIT,'(A,I4,A)') 'Taboo lists contain the lowest ',NTAB,' minima'
920:       ELSE IF (TABOOT) THEN936:       ELSE IF (TABOOT) THEN
921:          WRITE(MYUNIT,'(A,I4,A)') 'Taboo list contains the lowest ',NTAB,' minima'937:          WRITE(MYUNIT,'(A,I4,A)') 'Taboo list contains the lowest ',NTAB,' minima'
922:       ENDIF938:       ENDIF


r33130/minpermdist.f90 2017-08-07 15:30:56.185598380 +0100 r33129/minpermdist.f90 2017-08-07 15:30:57.833620309 +0100
 76: DOUBLE PRECISION TEMPCOORDSA(DEGFREEDOMS), TEMPCOORDSB(DEGFREEDOMS) 76: DOUBLE PRECISION TEMPCOORDSA(DEGFREEDOMS), TEMPCOORDSB(DEGFREEDOMS)
 77: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES) 77: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES)
 78: SAVE NORBIT1, NORBIT2 78: SAVE NORBIT1, NORBIT2
 79: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS) 79: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS)
 80: CHARACTER(LEN=5) ZSYMSAVE 80: CHARACTER(LEN=5) ZSYMSAVE
 81: COMMON /SYS/ ZSYMSAVE 81: COMMON /SYS/ ZSYMSAVE
 82:  82: 
 83: DOUBLE PRECISION :: DINV 83: DOUBLE PRECISION :: DINV
 84:  84: 
 85: DISTANCE = 0.0D0 85: DISTANCE = 0.0D0
  86: RMAT=0.0D0
  87: !cv320: commented out this write statement, looks like it was just for debugging.
  88: ! WRITE(*,*) 'RMAT', RMAT
  89: ! sf344 temporary debug
  90: !      CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGID,DEBUG,RMAT)
  91: !      RMATBEST = RMAT
  92: !   RETURN
  93: ! end sf344
 86:  94: 
 87: ! jwrm2 95: ! jwrm2
 88: IF (GTHOMSONT) THEN 96: IF (GTHOMSONT) THEN
 89:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,DISTANCE,DIST2,RMATBEST) 97:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,DISTANCE,DIST2,RMATBEST)
 90:    RETURN 98:    RETURN
 91: ENDIF 99: ENDIF
 92: 100: 
 93: IF (RIGIDINIT) THEN101: IF (RIGIDINIT) THEN
 94:     IF(DEBUG) THEN102:     IF(DEBUG) THEN
 95:         IF(.NOT.(ANY(ABS(COORDSA(DEGFREEDOMS+1:3*NATOMS)) .GT. 1.0E-10))) THEN103:         IF(.NOT.(ANY(ABS(COORDSA(DEGFREEDOMS+1:3*NATOMS)) .GT. 1.0E-10))) THEN
687:       WRITE(MYUNIT,'(2(A,F20.10))') 'minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &695:       WRITE(MYUNIT,'(2(A,F20.10))') 'minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &
688:   &                         ' should be ',SQRT(DISTANCE)696:   &                         ' should be ',SQRT(DISTANCE)
689:    ENDIF697:    ENDIF
690: 698: 
691:    IF ((NFREEZE.GT.0).OR.MKTRAPT) THEN699:    IF ((NFREEZE.GT.0).OR.MKTRAPT) THEN
692:       RMATBEST(1:3,1:3)=0.0D0700:       RMATBEST(1:3,1:3)=0.0D0
693:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0701:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
694:    ELSE702:    ELSE
695:       RMATBEST=MATMUL(ROTINVB,RMATBEST)703:       RMATBEST=MATMUL(ROTINVB,RMATBEST)
696:    ENDIF704:    ENDIF
697:    IF (DEBUG) THEN 705: !  IF (DEBUG) THEN 
698:       WRITE(MYUNIT, '(A)') 'RMATBEST:'706: !     WRITE(MYUNIT, '(A)') 'RMATBEST:'
699:       WRITE(MYUNIT, '(3F20.10)') RMATBEST(1:3,1:3)707: !     WRITE(MYUNIT, '(3F20.10)') RMATBEST(1:3,1:3)
700:    ENDIF708: !  ENDIF
701:    COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!709:    COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!
702: 710: 
703: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!711: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
704: !  DO J1=1,(NATOMS/2)712: !  DO J1=1,(NATOMS/2)
705: !     XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-COORDSA(3*(J1-1)+1))**2+ &713: !     XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-COORDSA(3*(J1-1)+1))**2+ &
706: ! &                 (COORDSB(3*(J1-1)+2)-COORDSA(3*(J1-1)+2))**2+ &714: ! &                 (COORDSB(3*(J1-1)+2)-COORDSA(3*(J1-1)+2))**2+ &
707: ! &                 (COORDSB(3*(J1-1)+3)-COORDSA(3*(J1-1)+3))**2715: ! &                 (COORDSB(3*(J1-1)+3)-COORDSA(3*(J1-1)+3))**2
708: !  ENDDO716: !  ENDDO
709: !  PRINT '(A,F20.10)','XDUMMY=',XDUMMY717: !  PRINT '(A,F20.10)','XDUMMY=',XDUMMY
710: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!718: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


r33130/modcudalbfgs.F90 2017-08-07 15:30:56.401601259 +0100 r33129/modcudalbfgs.F90 2017-08-07 15:30:58.053623235 +0100
 58:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration 58:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration
 59:  59: 
 60:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix 60:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix
 61:  61: 
 62:         REAL(KIND=C_DOUBLE), DIMENSION(N), INTENT(INOUT) :: C_XCOORDS ! Coordinates 62:         REAL(KIND=C_DOUBLE), DIMENSION(N), INTENT(INOUT) :: C_XCOORDS ! Coordinates
 63:  63: 
 64:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy 64:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy
 65:                                             C_RMS, & ! RMS force 65:                                             C_RMS, & ! RMS force
 66:                                             POTENTIALTIME ! Time taken in calculating potential - not used in GMIN 66:                                             POTENTIALTIME ! Time taken in calculating potential - not used in GMIN
 67:  67: 
 68:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info. 68:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info
 69:                                             C_CUDATIMET, & ! If true, print timing info.  69:                                             C_CUDATIMET, & ! If true, print timing info
 70:                                             C_ATOMRIGIDCOORDT, & ! If false, use rigid body coordinates 70:                                             C_ATOMRIGIDCOORDT, & ! If false, use rigid body coordinates
 71:                                             C_FREEZE, & ! If true, freeze some specified atoms 71:                                             C_FREEZE, & ! If true, freeze some specified atoms
 72:                                             PROJECT, & ! PROJECT is OPTIM only, always false 72:                                             PROJECT, & ! PROJECT is OPTIM only, always false
 73:                                             C_AACONVERGENCET ! If true, use more accurate method of calculating rigid RMS force 73:                                             C_AACONVERGENCET ! If true, use more accurate method of calculating rigid RMS force
 74:  74: 
 75:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3), INTENT(IN) :: C_FROZEN ! Logical array specifying frozen atoms 75:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3), INTENT(IN) :: C_FROZEN ! Logical array specifying frozen atoms
 76:  76: 
 77:         LOGICAL(KIND=C_BOOL), INTENT(OUT) :: C_MFLAG, & ! True if quench converged 77:         LOGICAL(KIND=C_BOOL), INTENT(OUT) :: C_MFLAG, & ! True if quench converged
 78:                                              C_COLDFUSION ! Set to true during minimization if cold fusion diagnosed 78:                                              C_COLDFUSION ! Set to true during minimization if cold fusion diagnosed
 79:  79: 
 80:         CHARACTER(LEN=1, KIND=C_CHAR), INTENT(IN) :: C_CUDAPOT ! Character specifying the CUDA potential to be used 80:         CHARACTER(LEN=1, KIND=C_CHAR), INTENT(IN) :: C_CUDAPOT ! Character specifying the CUDA potential to be used
 81:  81: 
 82:     END SUBROUTINE CUDA_LBFGS 82:     END SUBROUTINE CUDA_LBFGS
 83: END INTERFACE 83: END INTERFACE
 84:  84: 
 85: INTERFACE 85: INTERFACE
 86:     SUBROUTINE CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, &  86:     SUBROUTINE CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT, C_CUDAPOT, & 
 87:                                      C_LJADDATT) BIND(C,NAME="gminoptim_enegrad_cputogpu") 87:                                     C_CUDATIMET, POTENTIALTIME) BIND(C,NAME="setup_potential_cputogpu")
 88:  88: 
 89:         IMPORT :: C_INT, C_DOUBLE 89:         IMPORT :: C_INT, C_DOUBLE, C_BOOL, C_CHAR
 90:  90: 
 91:         INTEGER(KIND=C_INT), INTENT(IN) :: NATOMS, & ! No. of atoms 91:         INTEGER(KIND=C_INT), INTENT(IN) :: NATOMS, & ! No. of atoms
 92:                                            C_NADDTARGET ! Target cluster size (addressability) 92:                                            C_NADDTARGET ! Target cluster size (addressability)
 93:  93: 
 94:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(IN) :: COORDS ! Atomic coordinates 94:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(IN) :: COORDS ! Atomic coordinates
 95:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix 95:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix
 96:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_TOTENERGY ! Total energy of the system 96:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_TOTENERGY, & ! Total energy of the system
  97:                                             POTENTIALTIME ! Time taken in calculating potential - not used in GMIN
 97:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(OUT) :: C_GRADIENTS ! Gradient of the energy w.r.t. each atomic coordinate 98:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(OUT) :: C_GRADIENTS ! Gradient of the energy w.r.t. each atomic coordinate
 98:  99: 
 100:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_CUDATIMET ! If true, print timing info
 101: 
 102:         CHARACTER(LEN=1, KIND=C_CHAR), INTENT(IN) :: C_CUDAPOT ! Character specifying the CUDA potential to be used
 103: 
 99:     END SUBROUTINE CUDA_ENEGRAD_CPUTOGPU    104:     END SUBROUTINE CUDA_ENEGRAD_CPUTOGPU    
100: END INTERFACE105: END INTERFACE
101: #endif /* DUMMY_CUDA */106: #endif /* DUMMY_CUDA */
102: 107: 
103: CONTAINS108: CONTAINS
104: 109: 
105:     SUBROUTINE CUDA_LBFGS_WRAPPER(N, MUPDATE, XCOORDS, EPS, MFLAG, ENERGY, ITMAX, ITDONE, RESET)110:     SUBROUTINE CUDA_LBFGS_WRAPPER(N, MUPDATE, XCOORDS, EPS, MFLAG, ENERGY, ITMAX, ITDONE, RESET)
106: 111: 
107:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_LBFGS are converted directly112:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_LBFGS are converted directly
108: #ifndef DUMMY_CUDA113: #ifndef DUMMY_CUDA
272:         END IF277:         END IF
273: 278: 
274: #endif /* DUMMY_CUDA */279: #endif /* DUMMY_CUDA */
275: 280: 
276:     END SUBROUTINE CUDA_LBFGS_WRAPPER281:     END SUBROUTINE CUDA_LBFGS_WRAPPER
277: 282: 
278:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)283:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)
279: 284: 
280:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET285:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET
281: 286: 
282:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY287:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY, POTENTIALTIME
283:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS 
284: #ifndef DUMMY_CUDA288: #ifndef DUMMY_CUDA
 289:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS
285:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT290:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT
286: #else /* ifdef DUMMY_CUDA */291: #else /* ifdef DUMMY_CUDA */
 292:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: COORDS, C_GRADIENTS
287:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT293:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT
288: #endif294: #endif
 295: 
 296:         LOGICAL(KIND=C_BOOL) :: C_CUDATIMET
 297: 
 298:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
 299: 
289:         INTEGER :: X, I, J300:         INTEGER :: X, I, J
290: 301: 
291:         DOUBLE PRECISION :: TOTENERGY302:         DOUBLE PRECISION :: TOTENERGY
292:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS303:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS
293: 304: 
294: #ifndef DUMMY_CUDA305: #ifndef DUMMY_CUDA
295: 306: 
296:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN307:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN
297:             DO J = 1,NADDTARGET308:             DO J = 1,NADDTARGET
298:                 DO I = 1,NADDTARGET309:                 DO I = 1,NADDTARGET
299:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J)310:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J)
300:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J)311:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J)
301:                 END DO312:                 END DO
302:             END DO313:             END DO
303:         ELSE314:         ELSE
304:             C_LJADDREP(:) = 1.0315:             C_LJADDREP(:) = 1.0
305:             C_LJADDATT(:) = 1.0316:             C_LJADDATT(:) = 1.0
306:         END IF317:         END IF
307: 318: 
308:         C_NADDTARGET = NADDTARGET319:         C_NADDTARGET = NADDTARGET
 320:         C_CUDAPOT = CUDAPOT
 321:         C_CUDATIMET = CUDATIMET
309: 322: 
310:         ! Calculates the energy and gradients on the GPU using the GB potential323:         ! Calculates the energy and gradients on the GPU using the GB potential
311:         CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT)324:         CALL CUDA_ENEGRAD_CPUTOGPU(3*NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT, C_CUDAPOT, & 
 325:                                   C_CUDATIMET, POTENTIALTIME)
312: 326: 
313:         TOTENERGY = DBLE(C_TOTENERGY)327:         TOTENERGY = DBLE(C_TOTENERGY)
314:         328: 
315:         DO X = 1,(3*NATOMS)329:         DO X = 1,(3*NATOMS)
316:             GRADIENTS(X) = DBLE(C_GRADIENTS(X))330:             GRADIENTS(X) = DBLE(C_GRADIENTS(X))
317:         END DO331:         END DO
318: #endif332: #endif
319:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER333:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER
320: 334: 
321:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)335:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)
322:         IMPLICIT NONE336:         IMPLICIT NONE
323: 337: 
324:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET338:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET
325:         REAL(KIND=C_DOUBLE) :: C_ENERGY339:         REAL(KIND=C_DOUBLE) :: C_ENERGY, POTENTIALTIME
326:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)340:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)
327: #ifndef DUMMY_CUDA341: #ifndef DUMMY_CUDA
328:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT342:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT
329: #else /* ifdef DUMMY_CUDA */343: #else /* ifdef DUMMY_CUDA */
330:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT344:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT
331: #endif345: #endif
 346: 
 347:         LOGICAL(KIND=C_BOOL) :: C_CUDATIMET
 348: 
 349:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
 350: 
332:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)351:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)
333:         DOUBLE PRECISION    :: DELTA352:         DOUBLE PRECISION    :: DELTA
334:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)353:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)
335:         INTEGER             :: I, J354:         INTEGER             :: I, J
336: 355: 
337: #ifndef DUMMY_CUDA356: #ifndef DUMMY_CUDA
338:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN357:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN
339:             DO J = 1,NADDTARGET358:             DO J = 1,NADDTARGET
340:                 DO I = 1,NADDTARGET359:                 DO I = 1,NADDTARGET
341:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J)360:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J)
342:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J)361:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J)
343:                 END DO362:                 END DO
344:             END DO363:             END DO
345:         ELSE364:         ELSE
346:             C_LJADDREP(:) = 1.0365:             C_LJADDREP(:) = 1.0
347:             C_LJADDATT(:) = 1.0366:             C_LJADDATT(:) = 1.0
348:         END IF367:         END IF
349: 368: 
350:         C_NADDTARGET = NADDTARGET369:         C_NADDTARGET = NADDTARGET
 370:         C_CUDAPOT = CUDAPOT
 371:         C_CUDATIMET = CUDATIMET
351: 372: 
352:         DO I = 1, 3*NATOMS373:         DO I = 1, 3*NATOMS
353:             ! Plus374:             ! Plus
354:             COORDS(I) = COORDS(I) + DELTA375:             COORDS(I) = COORDS(I) + DELTA
355:             CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_ENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT)376:             CALL CUDA_ENEGRAD_CPUTOGPU(3*NATOMS, COORDS, C_ENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT, C_CUDAPOT, & 
 377:                                       C_CUDATIMET, POTENTIALTIME)
356:             GRAD_PLUS(:) = DBLE(C_GRADIENTS(:))378:             GRAD_PLUS(:) = DBLE(C_GRADIENTS(:))
357:             ! Minus379:             ! Minus
358:             COORDS(I) = COORDS(I) - 2.0D0 * DELTA380:             COORDS(I) = COORDS(I) - 2.0D0 * DELTA
359:             CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_ENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT)381:             CALL CUDA_ENEGRAD_CPUTOGPU(3*NATOMS, COORDS, C_ENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT, C_CUDAPOT, & 
 382:                                       C_CUDATIMET, POTENTIALTIME)
360:             GRAD_MINUS(:) = DBLE(C_GRADIENTS(:))383:             GRAD_MINUS(:) = DBLE(C_GRADIENTS(:))
361:             ! Reset coords384:             ! Reset coords
362:             COORDS(I) = COORDS(I) + DELTA385:             COORDS(I) = COORDS(I) + DELTA
363:             ! Calculate hessian386:             ! Calculate hessian
364:             HESSIAN(I, :) = (GRAD_PLUS(:) - GRAD_MINUS(:)) / (2.0D0 * DELTA)387:             HESSIAN(I, :) = (GRAD_PLUS(:) - GRAD_MINUS(:)) / (2.0D0 * DELTA)
365:         END DO388:         END DO
366: #endif389: #endif
367:     END SUBROUTINE CUDA_NUMERICAL_HESS390:     END SUBROUTINE CUDA_NUMERICAL_HESS
368: 391: 
369: END MODULE MODCUDALBFGS392: END MODULE MODCUDALBFGS


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0