hdiff output

r31980/amber12_interface.F90 2017-02-27 18:30:19.688090457 +0000 r31979/amber12_interface.F90 2017-02-27 18:30:21.892120092 +0000
 99:   end type pot_ene_rec_c 99:   end type pot_ene_rec_c
100: 100: 
101: !******************************************************************************101: !******************************************************************************
102: 102: 
103:   double precision, allocatable        :: amber12_masses(:)103:   double precision, allocatable        :: amber12_masses(:)
104:   type(amber12_atom), allocatable      :: amber12_atoms(:)104:   type(amber12_atom), allocatable      :: amber12_atoms(:)
105:   type(amber12_residue), allocatable   :: amber12_residues(:)105:   type(amber12_residue), allocatable   :: amber12_residues(:)
106: 106: 
107: contains107: contains
108: 108: 
109: subroutine amber12_setup(natoms, inpcrd_name_input, inpcrd_name_length) bind(C, name='amber12_setup')109: subroutine amber12_setup(natoms) bind(C, name='amber12_setup')
110: !110: !
111: ! Initialises all the stuff needed for AMBER to run... 111: ! Initialises all the stuff needed for AMBER to run... 
112: !112: !
113: ! Arguments113: ! Arguments
114: ! ---------114: ! ---------
115: !115: !
116: ! natoms(out): number of atoms from topology file116: ! natoms(out): number of atoms from topology file
117: ! inpcrd_name_input(in): file name for the inpcrd file. 
118: ! inpcrd_name_length(in): length of the file name 
119: !117: !
120:   use iso_c_binding, only: c_int, c_char, c_null_char118:   use iso_c_binding, only: c_int
121: #ifndef DUMMY_AMBER12119: #ifndef DUMMY_AMBER12
122:   use gb_alltasks_setup_mod120:   use gb_alltasks_setup_mod
123:   use pme_alltasks_setup_mod121:   use pme_alltasks_setup_mod
124: #ifdef DIRFRC_EFS122: #ifdef DIRFRC_EFS
125:   use ene_frc_splines_mod123:   use ene_frc_splines_mod
126: #endif /* DIRFRC_EFS */124: #endif /* DIRFRC_EFS */
127:   use file_io_dat_mod125:   use file_io_dat_mod
128:   use gb_ene_mod126:   use gb_ene_mod
129:   use inpcrd_dat_mod127:   use inpcrd_dat_mod
130:   use master_setup_mod128:   use master_setup_mod
142:   use remd_mod,       only : remd_method, bcast_remd_method, remd_setup, &140:   use remd_mod,       only : remd_method, bcast_remd_method, remd_setup, &
143:                              slave_remd_setup141:                              slave_remd_setup
144:   use remd_exchg_mod, only : setup_remd_randgen142:   use remd_exchg_mod, only : setup_remd_randgen
145:   use multipmemd_mod, only : setup_groups143:   use multipmemd_mod, only : setup_groups
146:   use pmemd_lib_mod144:   use pmemd_lib_mod
147:   use nmr_calls_mod145:   use nmr_calls_mod
148: #endif146: #endif
149: #endif /* DUMMY_AMBER12 */147: #endif /* DUMMY_AMBER12 */
150:   implicit none148:   implicit none
151: ! Arguments149: ! Arguments
152:   integer(c_int), intent(out)          :: natoms150:   integer(c_int), intent(out) :: natoms
153:   integer(c_int), intent(in)           :: inpcrd_name_length 
154:   character(kind=c_char, len=1), & 
155:             dimension(inpcrd_name_length), & 
156:             intent(in)                 :: inpcrd_name_input   
157: ! Variables  151: ! Variables  
158: #ifndef DUMMY_AMBER12152: #ifndef DUMMY_AMBER12
159:   double precision            :: max_erfc_relerr153:   double precision            :: max_erfc_relerr
160:   integer                     :: new_stack_limit    ! new stack limit154:   integer                     :: new_stack_limit    ! new stack limit
161:   integer                     :: num_ints = 0155:   integer                     :: num_ints = 0
162:   integer                     :: num_reals = 0156:   integer                     :: num_reals = 0
163:   integer                     :: i 
164:   logical                     :: terminal_flag = .false.157:   logical                     :: terminal_flag = .false.
165:   character(len=80)           :: inpcrd_name_string 
166: 158: 
167: #ifdef AMBMPI159: #ifdef AMBMPI
168: ! Establish pmemd communicators and process ranks160: ! Establish pmemd communicators and process ranks
169:   call setup_groups161:   call setup_groups
170: 162: 
171:   master = mytaskid .eq. 0 ! Make task 0 the master:163:   master = mytaskid .eq. 0 ! Make task 0 the master:
172: 164: 
173: ! CUDA version with GB can run on 1 GPU / task.165: ! CUDA version with GB can run on 1 GPU / task.
174: ! All other combinations require .gt. 1166: ! All other combinations require .gt. 1
175: #ifndef CUDA167: #ifndef CUDA
203:   inpcrd_name = 'coords.inpcrd'195:   inpcrd_name = 'coords.inpcrd'
204:   prmtop_name = 'coords.prmtop'196:   prmtop_name = 'coords.prmtop'
205:   restrt_name = 'min.rst' 197:   restrt_name = 'min.rst' 
206:   refc_name   = 'refc'198:   refc_name   = 'refc'
207:   mdvel_name  = 'mdvel'199:   mdvel_name  = 'mdvel'
208:   mden_name   = 'mden'200:   mden_name   = 'mden'
209:   mdcrd_name  = 'mdcrd'201:   mdcrd_name  = 'mdcrd'
210:   mdinfo_name = 'mdinfo'202:   mdinfo_name = 'mdinfo'
211:   logfile_name = 'amberlogfile'203:   logfile_name = 'amberlogfile'
212: 204: 
213: ! First check whether we're dealing with a C or Fortran string and set the 
214: ! length accordingly. 
215:   if (inpcrd_name_input(inpcrd_name_length) == c_null_char) then 
216:     do i = 1, inpcrd_name_length - 1 
217:       inpcrd_name_string(i:i) = inpcrd_name_input(i) 
218:     end do 
219:     inpcrd_name_string(inpcrd_name_length:80)     = " " 
220:   else 
221:     do i = 1, inpcrd_name_length 
222:       inpcrd_name_string(i:i) = inpcrd_name_input(i) 
223:     end do 
224:     inpcrd_name_string(inpcrd_name_length + 1:80) = " " 
225:   end if 
226:   inpcrd_name = inpcrd_name_string 
227:  
228:   if (master) then205:   if (master) then
229:     call master_setup(num_ints, num_reals, new_stack_limit, terminal_flag)206:     call master_setup(num_ints, num_reals, new_stack_limit, terminal_flag)
230: #ifdef CUDA207: #ifdef CUDA
231:   else208:   else
232:     !RCW: Call to set device is mostly redundant now due to209:     !RCW: Call to set device is mostly redundant now due to
233:     !     removal of -gpu command line argument. But leave for now.210:     !     removal of -gpu command line argument. But leave for now.
234:     call gpu_set_device(-1)211:     call gpu_set_device(-1)
235:     call gpu_init()212:     call gpu_init()
236: #ifdef AMBMPI213: #ifdef AMBMPI
237:     call gpu_send_slave_device_info()214:     call gpu_send_slave_device_info()


r31980/CMakeLists.txt 2017-02-27 18:30:19.908093415 +0000 r31979/CMakeLists.txt 2017-02-27 18:30:22.276125256 +0000
170: set_module_dir(extralib)170: set_module_dir(extralib)
171: set_target_properties(extralib PROPERTIES LINKER_LANGUAGE "Fortran") 171: set_target_properties(extralib PROPERTIES LINKER_LANGUAGE "Fortran") 
172: set_target_properties(extralib PROPERTIES COMPILE_DEFINITIONS "${COMPILE_DEFINITIONS};${DUMMY_CPP_FLAGS};__SPARSE")172: set_target_properties(extralib PROPERTIES COMPILE_DEFINITIONS "${COMPILE_DEFINITIONS};${DUMMY_CPP_FLAGS};__SPARSE")
173: 173: 
174: # Link CHOLMOD, if we're using it174: # Link CHOLMOD, if we're using it
175: option(WITH_SPARSE "Compile with SuiteSparse for sparse hessian diagonalisation" OFF)175: option(WITH_SPARSE "Compile with SuiteSparse for sparse hessian diagonalisation" OFF)
176: if(WITH_SPARSE)176: if(WITH_SPARSE)
177:    add_subdirectory(${SVN_ROOT}/SuiteSparse SuiteSparse)177:    add_subdirectory(${SVN_ROOT}/SuiteSparse SuiteSparse)
178:    include_directories(${SVN_ROOT}/SuiteSparse/SuiteSparse/CHOLMOD/Include)178:    include_directories(${SVN_ROOT}/SuiteSparse/SuiteSparse/CHOLMOD/Include)
179:    include_directories(${SVN_ROOT}/SuiteSparse/SuiteSparse/SuiteSparse_config)179:    include_directories(${SVN_ROOT}/SuiteSparse/SuiteSparse/SuiteSparse_config)
180:    if(SuiteSparse_GPU) 
181:       include_directories(${CUDA_TOOLKIT_ROOT_DIR}/include) 
182:    endif(SuiteSparse_GPU) 
183:    add_definitions(-D__SPARSE)180:    add_definitions(-D__SPARSE)
184:    file(GLOB SPARSE_INTERFACE sparse/*.f90181:    file(GLOB SPARSE_INTERFACE sparse/*.f90
185:                               sparse/*.c ) 182:                               sparse/*.c ) 
186:    list(APPEND GMIN_LIB_SOURCES ${SPARSE_INTERFACE})183:    list(APPEND GMIN_LIB_SOURCES ${SPARSE_INTERFACE})
187: endif(WITH_SPARSE)184: endif(WITH_SPARSE)
188: 185: 
189: # Make a gmin library186: # Make a gmin library
190: add_library(gminlib ${GMIN_LIB_SOURCES})187: add_library(gminlib ${GMIN_LIB_SOURCES})
191: set_module_dir(gminlib)188: set_module_dir(gminlib)
192: set_module_depends(gminlib extralib)189: set_module_depends(gminlib extralib)


r31980/commons.f90 2017-02-27 18:30:20.128096373 +0000 r31979/commons.f90 2017-02-27 18:30:22.540128806 +0000
597:       LOGICAL :: SELECTMOVET597:       LOGICAL :: SELECTMOVET
598:       DOUBLE PRECISION, ALLOCATABLE :: SELTRANSSTEP(:), SELROTSCALE(:), SELMOVPROB(:)598:       DOUBLE PRECISION, ALLOCATABLE :: SELTRANSSTEP(:), SELROTSCALE(:), SELMOVPROB(:)
599:       INTEGER, ALLOCATABLE :: SELMOVFREQ(:), SELBEGIN(:), SELEND(:), SELSIZE(:)599:       INTEGER, ALLOCATABLE :: SELMOVFREQ(:), SELBEGIN(:), SELEND(:), SELSIZE(:)
600:       INTEGER :: SELMOVNO600:       INTEGER :: SELMOVNO
601: 601: 
602: !QUIP variables602: !QUIP variables
603:       CHARACTER(LEN=3) :: QUIPATOMTYPE603:       CHARACTER(LEN=3) :: QUIPATOMTYPE
604:       CHARACTER(LEN=10240) :: QUIPARGSTR604:       CHARACTER(LEN=10240) :: QUIPARGSTR
605: 605: 
606: ! FEBH variables606: ! FEBH variables
607:       LOGICAL           :: FEBHT, SPARSET, SPARSE_BENCH607:       LOGICAL           :: FEBHT, SPARSET
608:       DOUBLE PRECISION  :: FEBH_POT_ENE, FETEMP, ZERO_THRESH608:       DOUBLE PRECISION  :: FEBH_POT_ENE, FETEMP, ZERO_THRESH
609:       DOUBLE PRECISION, PARAMETER :: SMALL_DOUBLE = 1.0D-100609:       DOUBLE PRECISION, PARAMETER :: SMALL_DOUBLE = 1.0D-100
610:       DOUBLE PRECISION, ALLOCATABLE :: QENERGIES(:), QCOORDINATES(:,:), QPE(:)610:       DOUBLE PRECISION, ALLOCATABLE :: QENERGIES(:), QCOORDINATES(:,:), QPE(:)
611:       INTEGER           :: FE_FILE_UNIT611:       INTEGER           :: FE_FILE_UNIT
612: ! Variables for converging to a certain separation of zero and non-zero eigenvalues.612: ! Variables for converging to a certain separation of zero and non-zero eigenvalues.
613:       DOUBLE PRECISION  :: MIN_ZERO_SEP613:       DOUBLE PRECISION  :: MIN_ZERO_SEP
614:       INTEGER           :: MAX_ATTEMPTS614:       INTEGER           :: MAX_ATTEMPTS
615: 615: 
616: !Character string for CUDA potential616: !Character string for CUDA potential
617:       CHARACTER(LEN=1) :: CUDAPOT617:       CHARACTER(LEN=1) :: CUDAPOT


r31980/countatoms.f90 2017-02-27 18:30:20.348099331 +0000 r31979/countatoms.f90 2017-02-27 18:30:22.760131764 +0000
 35:     INTEGER, INTENT(IN) :: NPAR 35:     INTEGER, INTENT(IN) :: NPAR
 36:     INTEGER :: EOF,NRES,SEQ(500),I_RES,NOGLY,GLY, MYUNIT,J 36:     INTEGER :: EOF,NRES,SEQ(500),I_RES,NOGLY,GLY, MYUNIT,J
 37:     LOGICAL :: YESNO, YESNOA, YESNOC, YESNOAMH, YESNOA9, YESNOOPEP 37:     LOGICAL :: YESNO, YESNOA, YESNOC, YESNOAMH, YESNOA9, YESNOOPEP
 38:     CHARACTER(LEN=5) TARFL 38:     CHARACTER(LEN=5) TARFL
 39:     CHARACTER(LEN=10)  CHECK 39:     CHARACTER(LEN=10)  CHECK
 40:     CHARACTER(LEN=80) MYLINE,INPCRD1 40:     CHARACTER(LEN=80) MYLINE,INPCRD1
 41:     LOGICAL GCBHT, END, GBHT 41:     LOGICAL GCBHT, END, GBHT
 42:     DOUBLE PRECISION GCMU 42:     DOUBLE PRECISION GCMU
 43:     INTEGER GCNATOMS, LUNIT, GETUNIT 43:     INTEGER GCNATOMS, LUNIT, GETUNIT
 44:     CHARACTER(LEN=16) WORD 44:     CHARACTER(LEN=16) WORD
 45:     CHARACTER(LEN=20) OSTRING 
 46:      45:     
 47: ! ds656> First scan 'data' for GCBH and GBH, since the latter  46: ! ds656> First scan 'data' for GCBH and GBH, since the latter 
 48: !        affects atom counting. 47: !        affects atom counting.
 49:  48: 
 50:     LUNIT=GETUNIT() 49:     LUNIT=GETUNIT()
 51:     OPEN(UNIT=LUNIT,FILE='data',STATUS='OLD') 50:     OPEN(UNIT=LUNIT,FILE='data',STATUS='OLD')
 52:      51:     
 53: 190 CALL INPUT(END,LUNIT) 52: 190 CALL INPUT(END,LUNIT)
 54:     IF (.NOT. END) THEN 53:     IF (.NOT. END) THEN
 55:        CALL READU(WORD) 54:        CALL READU(WORD)
145:           if (seq(i_res).eq.8) GLY = GLY +1144:           if (seq(i_res).eq.8) GLY = GLY +1
146:        enddo145:        enddo
147:        146:        
148:        Number_of_Atoms = NOGLY*3 + GLY*2147:        Number_of_Atoms = NOGLY*3 + GLY*2
149:     ELSE IF (YESNOA9) THEN148:     ELSE IF (YESNOA9) THEN
150:        ! OPEN(UNIT=7,FILE='coords.gayberne',STATUS='OLD')149:        ! OPEN(UNIT=7,FILE='coords.gayberne',STATUS='OLD')
151:        ! PRINT '(A)','reading coordinates from file coords.gayberne'150:        ! PRINT '(A)','reading coordinates from file coords.gayberne'
152:        ! inpcrd1='coords.inpcrd'151:        ! inpcrd1='coords.inpcrd'
153:        ! inpcrd1=trim(adjustl(inpcrd1))152:        ! inpcrd1=trim(adjustl(inpcrd1))
154:        ! call amberinterface(Number_of_Atoms,1,inpcrd1,MYUNIT)153:        ! call amberinterface(Number_of_Atoms,1,inpcrd1,MYUNIT)
155:        WRITE(OSTRING,'(A)') 'coords.inpcrd'154:        CALL AMBER12_SETUP(NUMBER_OF_ATOMS)
156:        CALL AMBER12_SETUP(NUMBER_OF_ATOMS, OSTRING, LEN(OSTRING)) 
157:        IF (NUMBER_OF_ATOMS .EQ. 0) THEN155:        IF (NUMBER_OF_ATOMS .EQ. 0) THEN
158:           NUMBER_OF_ATOMS = COUNT_ATOMS_FROM_INPCRD('coords.inpcrd')156:           NUMBER_OF_ATOMS = COUNT_ATOMS_FROM_INPCRD('coords.inpcrd')
159:        END IF157:        END IF
160:        158:        
161:     ELSEIF (YESNOC) THEN159:     ELSEIF (YESNOC) THEN
162:        OPEN(UNIT=7,FILE='input.crd',STATUS='OLD')160:        OPEN(UNIT=7,FILE='input.crd',STATUS='OLD')
163:        do161:        do
164:           read(7,*) myline162:           read(7,*) myline
165:           if (myline(1:1)=='*') then ! SAT This is the goddamn CHARMM comment line163:           if (myline(1:1)=='*') then ! SAT This is the goddamn CHARMM comment line
166:              cycle164:              cycle


r31980/getparams.f 2017-02-27 18:30:21.672117134 +0000 r31979/getparams.f 2017-02-27 18:30:24.184150911 +0000
 35:       COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR, 35:       COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR,
 36:      &                NERROR, IR, ECHO, LAST, CAT 36:      &                NERROR, IR, ECHO, LAST, CAT
 37:       DOUBLE PRECISION DUMMY 37:       DOUBLE PRECISION DUMMY
 38:       CHARACTER ZDUM*5 38:       CHARACTER ZDUM*5
 39:       INTEGER J2, NELEMENTS, LSYS, NTYPE(105), IOS, NARG, FILTH, FILTH2, NTYPES 39:       INTEGER J2, NELEMENTS, LSYS, NTYPE(105), IOS, NARG, FILTH, FILTH2, NTYPES
 40:       INTEGER, ALLOCATABLE :: ELEMENT_NUMS(:) 40:       INTEGER, ALLOCATABLE :: ELEMENT_NUMS(:)
 41:       CHARACTER FNAME*80, TSTRING*80 41:       CHARACTER FNAME*80, TSTRING*80
 42:       CHARACTER(LEN=80) :: SYS 42:       CHARACTER(LEN=80) :: SYS
 43:       CHARACTER WORD*16 43:       CHARACTER WORD*16
 44:       CHARACTER(LEN=10)  check 44:       CHARACTER(LEN=10)  check
 45:       CHARACTER(LEN=20) OTEMP, OSTRING, CSTRING, FIRSTAMBERSTR 45:       CHARACTER(LEN=20) OTEMP, OSTRING, CSTRING
 46:       CHARACTER(LEN=21) DSTRING1, DSTRING2 46:       CHARACTER(LEN=21) DSTRING1, DSTRING2
 47:       CHARACTER(LEN=80) ARGSTRING, MYLINE 47:       CHARACTER(LEN=80) ARGSTRING, MYLINE
 48:       LOGICAL AMH 48:       LOGICAL AMH
 49:       INTEGER :: NRES,I_RES,NOGLY 49:       INTEGER :: NRES,I_RES,NOGLY
 50:       INTEGER LUNIT,GETUNIT, EPOCH 50:       INTEGER LUNIT,GETUNIT, EPOCH
 51:  51: 
 52: ! TIME LIMITING BINARIES - uncomment the below 52: ! TIME LIMITING BINARIES - uncomment the below
 53: ! csw34> Workshop time limiting 53: ! csw34> Workshop time limiting
 54: ! find the number of seconds since 0000 on 1/1/1970 (unix epoch) 54: ! find the number of seconds since 0000 on 1/1/1970 (unix epoch)
 55: !      EPOCH=TIME() 55: !      EPOCH=TIME()
255:          ENDDO255:          ENDDO
256: 211      CONTINUE256: 211      CONTINUE
257:       ELSE IF (WORD.EQ.'RINGPOLYMER') THEN257:       ELSE IF (WORD.EQ.'RINGPOLYMER') THEN
258:          RINGPOLYMERT=.TRUE.258:          RINGPOLYMERT=.TRUE.
259:          GOTO 200259:          GOTO 200
260:       ELSE IF (WORD.EQ.'VARIABLES') THEN260:       ELSE IF (WORD.EQ.'VARIABLES') THEN
261:          VARIABLES=.TRUE.261:          VARIABLES=.TRUE.
262:          GOTO 200262:          GOTO 200
263:       ELSE IF (WORD .EQ. 'AMBER12') THEN263:       ELSE IF (WORD .EQ. 'AMBER12') THEN
264:          AMBER12T = .TRUE.264:          AMBER12T = .TRUE.
265:          IF (NITEMS == 3) THEN265:          CALL AMBER12_SETUP(NATOMS)
266:             CALL READA(FIRSTAMBERSTR) 
267:             IF (FILTH2 .NE. 0) THEN 
268:                WRITE(OTEMP, *) FILTH2 
269:                WRITE(OSTRING,'(A)') TRIM(ADJUSTL(FIRSTAMBERSTR))//'.'//TRIM(ADJUSTL(OTEMP)) 
270:                WRITE(*,*) 'ostring=', OSTRING 
271:             ELSE 
272:                WRITE(OSTRING,'(A)') TRIM(ADJUSTL(FIRSTAMBERSTR)) 
273:             END IF 
274:             CALL READA(FIRSTAMBERSTR) 
275:             IF (TRIM(ADJUSTL(FIRSTAMBERSTR)) /= 'inpcrd') THEN 
276:                WRITE(OSTRING,'(A)') 'coords.inpcrd' 
277:             END IF 
278:          END IF 
279:          CALL AMBER12_SETUP(NATOMS, OSTRING, LEN(OSTRING)) 
280: ! Go to the end of the routine once we have the correct number of atoms.266: ! Go to the end of the routine once we have the correct number of atoms.
281:          GOTO 400267:          GOTO 400
282:       ELSE IF (WORD.EQ.'AMBER9'.OR.(WORD.EQ.'NAB')) THEN268:       ELSE IF (WORD.EQ.'AMBER9'.OR.(WORD.EQ.'NAB')) THEN
283:          IF(WORD.EQ.'AMBER9') AMBERT=.TRUE.269:          IF(WORD.EQ.'AMBER9') AMBERT=.TRUE.
284:          IF(WORD.EQ.'NAB') NABT=.TRUE.270:          IF(WORD.EQ.'NAB') NABT=.TRUE.
285: !         NATOMS=0271: !         NATOMS=0
286:         INPCRD='coords.inpcrd'272:         INPCRD='coords.inpcrd'
287:        IF(NITEMS<3) then273:        IF(NITEMS<3) then
288:          CALL READA(amberstr)274:          CALL READA(amberstr)
289:           IF (FILTH2.NE.0) THEN275:           IF (FILTH2.NE.0) THEN


r31980/keywords.f 2017-02-27 18:30:20.576102397 +0000 r31979/keywords.f 2017-02-27 18:30:22.988134829 +0000
1114:       FEBHT = .FALSE.1114:       FEBHT = .FALSE.
1115:       FETEMP = 0.0D01115:       FETEMP = 0.0D0
1116: ! khs26> This requires a minimum separation between the zero and non-zero1116: ! khs26> This requires a minimum separation between the zero and non-zero
1117: ! eigenvalues for runs using free energy basin-hopping. This is based on1117: ! eigenvalues for runs using free energy basin-hopping. This is based on
1118: ! the minimum found in quench zero.1118: ! the minimum found in quench zero.
1119:       MIN_ZERO_SEP = 0.0D01119:       MIN_ZERO_SEP = 0.0D0
1120:       MAX_ATTEMPTS = 51120:       MAX_ATTEMPTS = 5
1121: ! Use sparse hessian methods1121: ! Use sparse hessian methods
1122:       SPARSET = .FALSE.1122:       SPARSET = .FALSE.
1123:       ZERO_THRESH = 0.0D01123:       ZERO_THRESH = 0.0D0
1124:       SPARSE_BENCH = .FALSE. 
1125: ! khs26> Rotamer move stuff1124: ! khs26> Rotamer move stuff
1126:       ROTAMER_MOVET = .FALSE.1125:       ROTAMER_MOVET = .FALSE.
1127: 1126: 
1128: 1127: 
1129:       CUDAT=.FALSE.1128:       CUDAT=.FALSE.
1130:       CUDAPOT=' '1129:       CUDAPOT=' '
1131:       CUDATIMET=.FALSE.1130:       CUDATIMET=.FALSE.
1132:       GCBHT=.FALSE.1131:       GCBHT=.FALSE.
1133:       SEMIGRAND_MUT=.FALSE.1132:       SEMIGRAND_MUT=.FALSE.
1134:       USEROT=.FALSE.1133:       USEROT=.FALSE.
1991:       ELSE IF (WORD.EQ.'CG') THEN1990:       ELSE IF (WORD.EQ.'CG') THEN
1992:          LBFGST=.FALSE.1991:          LBFGST=.FALSE.
1993:          CONJG=.TRUE.1992:          CONJG=.TRUE.
1994: !1993: !
1995: ! sf344> Start of AMBER-related keywords1994: ! sf344> Start of AMBER-related keywords
1996: !1995: !
1997:       ELSE IF (WORD.EQ.'AMBERMDSTEPS') THEN1996:       ELSE IF (WORD.EQ.'AMBERMDSTEPS') THEN
1998:         MDSTEPT = .TRUE.1997:         MDSTEPT = .TRUE.
1999:       ELSE IF (WORD .EQ. 'AMBER12') THEN1998:       ELSE IF (WORD .EQ. 'AMBER12') THEN
2000:         AMBER12T = .TRUE.1999:         AMBER12T = .TRUE.
2001:         SETCHIRAL = .TRUE. 
2002:         IF(.NOT.ALLOCATED(COORDS1)) ALLOCATE(COORDS1(3*NATOMS))2000:         IF(.NOT.ALLOCATED(COORDS1)) ALLOCATE(COORDS1(3*NATOMS))
2003:         IF(ALLOCATED(COORDS)) DEALLOCATE(COORDS)2001:         IF(ALLOCATED(COORDS)) DEALLOCATE(COORDS)
2004:         ! Read the coords from AMBER12 into COORDS1(:)2002:         ! Read the coords from AMBER12 into COORDS1(:)
2005:         CALL AMBER12_GET_COORDS(NATOMS, COORDS1(:))2003:         CALL AMBER12_GET_COORDS(NATOMS, COORDS1(:))
2006:         ALLOCATE(COORDS(3*NATOMS,NPAR))2004:         ALLOCATE(COORDS(3*NATOMS,NPAR))
2007:         DO J1=1,NPAR2005:         DO J1=1,NPAR
2008:            COORDS(:,J1) = COORDS1(:)2006:            COORDS(:,J1) = COORDS1(:)
2009:         END DO2007:         END DO
2010: 2008: 
2011:       ELSE IF (WORD.EQ.'AMBER9') THEN2009:       ELSE IF (WORD.EQ.'AMBER9') THEN
6058:          CALL READI(SOFT_SPHERE_NTYPEA)6056:          CALL READI(SOFT_SPHERE_NTYPEA)
6059:          NTYPEA = SOFT_SPHERE_NTYPEA6057:          NTYPEA = SOFT_SPHERE_NTYPEA
6060: 6058: 
6061:       ELSE IF (WORD.EQ.'SORT') THEN6059:       ELSE IF (WORD.EQ.'SORT') THEN
6062:          SORTT=.TRUE.6060:          SORTT=.TRUE.
6063: 6061: 
6064:       ELSE IF (WORD.EQ.'SPARSE') THEN6062:       ELSE IF (WORD.EQ.'SPARSE') THEN
6065:          SPARSET=.TRUE.6063:          SPARSET=.TRUE.
6066:          IF (NITEMS.GT.1) CALL READF(ZERO_THRESH)6064:          IF (NITEMS.GT.1) CALL READF(ZERO_THRESH)
6067: 6065: 
6068:       ELSE IF (WORD .EQ. 'SPARSEBENCH') THEN 
6069:          SPARSE_BENCH=.TRUE. 
6070:  
6071: !     ELSE IF (WORD.EQ.'SQUEEZE') THEN6066: !     ELSE IF (WORD.EQ.'SQUEEZE') THEN
6072: !        CALL READI(NVEC)6067: !        CALL READI(NVEC)
6073: !        SQUEEZET=.TRUE.6068: !        SQUEEZET=.TRUE.
6074: !        IF (NITEMS.GT.2) THEN6069: !        IF (NITEMS.GT.2) THEN
6075: !           CALL READF(XX)6070: !           CALL READF(XX)
6076: !           SQUEEZER=XX6071: !           SQUEEZER=XX
6077: !        ENDIF6072: !        ENDIF
6078: !        IF (NITEMS.GT.3) THEN6073: !        IF (NITEMS.GT.3) THEN
6079: !           CALL READF(XX)6074: !           CALL READF(XX)
6080: !           SQUEEZED=XX6075: !           SQUEEZED=XX


r31980/potential.f90 2017-02-27 18:30:20.792105301 +0000 r31979/potential.f90 2017-02-27 18:30:23.304139078 +0000
   1: 
  1: !   GMIN: A program for finding global minima  2: !   GMIN: A program for finding global minima
  2: !   Copyright (C) 1999-2006 David J. Wales  3: !   Copyright (C) 1999-2006 David J. Wales
  3: !   This file is part of GMIN.  4: !   This file is part of GMIN.
  4: !  5: !
  5: !   GMIN is free software; you can redistribute it and/or modify  6: !   GMIN is free software; you can redistribute it and/or modify
  6: !   it under the terms of the GNU General Public License as published by  7: !   it under the terms of the GNU General Public License as published by
  7: !   the Free Software Foundation; either version 2 of the License, or  8: !   the Free Software Foundation; either version 2 of the License, or
  8: !   (at your option) any later version.  9: !   (at your option) any later version.
  9: ! 10: !
 10: !   GMIN is distributed in the hope that it will be useful, 11: !   GMIN is distributed in the hope that it will be useful,
477:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)478:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)
478:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)479:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
479:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0480:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
480:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)481:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
481:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0482:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
482:       END IF483:       END IF
483: !copied from OPTIM, rbody part not tested484: !copied from OPTIM, rbody part not tested
484:       IF (SECT) THEN485:       IF (SECT) THEN
485:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)486:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
486:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))487:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))
487:          IF (CUDAT) THEN488:          CALL AMBER12_NUM_HESS(NATOMS,X, DELTA=1.0D-5, HESSIAN=HESS(:, :))
488:             CALL CUDA_NUMERICAL_HESS(NATOMS, X, HESS, DELTA=1.0D-6) 
489:          ELSE 
490:             CALL AMBER12_NUM_HESS(NATOMS,X, DELTA=1.0D-5, HESSIAN=HESS(:, :)) 
491:          END IF 
492:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT) ) THEN489:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT) ) THEN
493:             CALL TRANSFORMHESSIAN(HESS, GRADATOMS, XRIGIDCOORDS,XRIGIDHESS, RBAANORMALMODET)490:             CALL TRANSFORMHESSIAN(HESS, GRADATOMS, XRIGIDCOORDS,XRIGIDHESS, RBAANORMALMODET)
494:             HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0491:             HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0
495:             HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0492:             HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0
496:             HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)493:             HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)
497:          END IF494:          END IF
498:       END IF495:       END IF
499: ! AMBER 9 Energy and gradient calls496: ! AMBER 9 Energy and gradient calls
500:    ELSE IF (AMBERT) THEN497:    ELSE IF (AMBERT) THEN
501: ! hk286 > Generalised rigid body498: ! hk286 > Generalised rigid body
534: 531: 
535:       IF (SECT) THEN532:       IF (SECT) THEN
536:          VNEW=0.0D0533:          VNEW=0.0D0
537: ! khs26> Copied analytical second derivatives from OPTIM534: ! khs26> Copied analytical second derivatives from OPTIM
538:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN535:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
539:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)536:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
540:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)537:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
541:          END IF538:          END IF
542:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)539:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
543:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))540:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))
544:          IF (.NOT. ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))541:          IF (CUDAT) THEN
545:          TEMPHESS(:)=0.0D0542:             CALL CUDA_NUMERICAL_HESS(NATOMS, X, HESS, DELTA=1.0D-6)
546:          CALL MME2WRAPPER(X, ENERGY, VNEW, TEMPHESS, ATMASS1, GRAD1)543:          ELSE
547:          IF (SQRT(ABS(ENERGY)) > 1.0D40) THEN544:             IF (.NOT. ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))
548:             EREAL = ENERGY545:             TEMPHESS(:)=0.0D0
549:             WRITE(MYUNIT, *) "Linear dihedral detected in NAB routines (sff2.c)."546:             CALL MME2WRAPPER(X, ENERGY, VNEW, TEMPHESS, ATMASS1, GRAD1)
550:          END IF547:             IF (SQRT(ABS(ENERGY)) > 1.0D40) THEN
551:          K=1548:                EREAL = ENERGY
552:          DO I=1, 3*NATOMS549:                WRITE(MYUNIT, *) "Linear dihedral detected in NAB routines (sff2.c)."
553:             DO J=1, 3*NATOMS550:             END IF
554:                HESS(I, J)=TEMPHESS(K)551:             K=1
555:                K=K + 1552:             DO I=1, 3*NATOMS
 553:                DO J=1, 3*NATOMS
 554:                   HESS(I, J)=TEMPHESS(K)
 555:                   K=K + 1
 556:                END DO
556:             END DO557:             END DO
557:          END DO558:             DEALLOCATE(TEMPHESS)
558:          DEALLOCATE(TEMPHESS)559:          END IF
559:          IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN560:          IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
560:             CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)561:             CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)
561:             CALL TRANSFORMGRAD(VNEW, XRIGIDCOORDS, XRIGIDGRAD)562:             CALL TRANSFORMGRAD(VNEW, XRIGIDCOORDS, XRIGIDGRAD)
562:             X(DEGFREEDOMS+1:3*NATOMS)=0.0D0563:             X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
563:             X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)564:             X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
564:             VNEW(DEGFREEDOMS+1:3*NATOMS)=0.0D0565:             VNEW(DEGFREEDOMS+1:3*NATOMS)=0.0D0
565:             VNEW(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)566:             VNEW(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
566:             HESS(DEGFREEDOMS+1:3*NATOMS,:)=0.0D0567:             HESS(DEGFREEDOMS+1:3*NATOMS,:)=0.0D0
567:             HESS(:, DEGFREEDOMS+1:3*NATOMS)=0.0D0568:             HESS(:, DEGFREEDOMS+1:3*NATOMS)=0.0D0
568:             HESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)=XRIGIDHESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)569:             HESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)=XRIGIDHESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)


r31980/quench.F 2017-02-27 18:30:21.016108314 +0000 r31979/quench.F 2017-02-27 18:30:23.524142036 +0000
318:          ELSE318:          ELSE
319: ! Added timing for LBFGS call319: ! Added timing for LBFGS call
320:             CALL CPU_TIME(T_LBFGS_START)320:             CALL CPU_TIME(T_LBFGS_START)
321:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)321:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)
322:             CALL CPU_TIME(T_LBFGS_END)322:             CALL CPU_TIME(T_LBFGS_END)
323:             IF (FEBHT.AND.DEBUG) WRITE(MYUNIT, '(A, F10.3)') 'Time to minimise (s):', T_LBFGS_END - T_LBFGS_START323:             IF (FEBHT.AND.DEBUG) WRITE(MYUNIT, '(A, F10.3)') 'Time to minimise (s):', T_LBFGS_END - T_LBFGS_START
324:          ENDIF324:          ENDIF
325:          IF (EVAPREJECT) RETURN325:          IF (EVAPREJECT) RETURN
326:       ENDIF326:       ENDIF
327:       !327:       !
328:       IF (FEBHT .AND. CFLAG .AND. (.NOT. SPARSE_BENCH)) THEN328:       IF (FEBHT .AND. CFLAG) THEN
329:       ! Calculate the free energy329:       ! Calculate the free energy
330:          NUM_ZERO_EVS=6330:          NUM_ZERO_EVS=6
331:          IF (NATOMS.EQ.2) NUM_ZERO_EVS=5331:          IF (NATOMS.EQ.2) NUM_ZERO_EVS=5
332:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)332:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
333:          ALLOCATE(HESS(3*NATOMS,3*NATOMS))333:          ALLOCATE(HESS(3*NATOMS,3*NATOMS))
334:          CALL POTENTIAL(P,GRAD,EREAL,.TRUE.,.TRUE.)334:          CALL POTENTIAL(P,GRAD,EREAL,.TRUE.,.TRUE.)
335:          CALL MASSWT()335:          CALL MASSWT()
336: #ifdef __SPARSE336: #ifdef __SPARSE
337:          IF (SPARSET) THEN337:          IF (SPARSET) THEN
338:              IF (DEBUG) THEN338:              IF (DEBUG) THEN
475:             ITDET=LOG(ITDET)/2.0D0 475:             ITDET=LOG(ITDET)/2.0D0 
476:          ENDIF476:          ENDIF
477:          FEBH_POT_ENE=EREAL477:          FEBH_POT_ENE=EREAL
478: ! khs26> At zero temperature, free energy is potential energy. However, DLOG(0)=NaN478: ! khs26> At zero temperature, free energy is potential energy. However, DLOG(0)=NaN
479:          479:          
480:          IF (FETEMP .GT. 0.D0) THEN480:          IF (FETEMP .GT. 0.D0) THEN
481:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Potential energy=', FEBH_POT_ENE481:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Potential energy=', FEBH_POT_ENE
482:             EREAL=-FETEMP*(SYMFCTR+LOG(1.0D0/HORDER) - EREAL/FETEMP + (3*NATOMS-NUM_ZERO_EVS)*DLOG(FETEMP) - LOG_PROD_EV + ITDET)482:             EREAL=-FETEMP*(SYMFCTR+LOG(1.0D0/HORDER) - EREAL/FETEMP + (3*NATOMS-NUM_ZERO_EVS)*DLOG(FETEMP) - LOG_PROD_EV + ITDET)
483:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Harmonic superposition contribution=', EREAL - FEBH_POT_ENE483:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Harmonic superposition contribution=', EREAL - FEBH_POT_ENE
484:          END IF484:          END IF
485:       ELSE IF (FEBHT .AND. (.NOT. CFLAG) .AND. (.NOT. SPARSE_BENCH)) THEN485:       ELSE IF (FEBHT .AND. (.NOT. CFLAG)) THEN
486:          WRITE(MYUNIT, '(A)') 'Quench did not converge, not calculating free energy and adding 1E10 to energy.'486:          WRITE(MYUNIT, '(A)') 'Quench did not converge, not calculating free energy and adding 1E10 to energy.'
487:          EREAL=EREAL + 1.0D10487:          EREAL=EREAL + 1.0D10
488: ! Benchmarking of FEBH 
489:       ELSE IF (FEBHT .AND. CFLAG .AND. SPARSE_BENCH) THEN 
490:          CALL BENCH_SPARSE(P, NQ(NP), ZERO_THRESH, EREAL, LOG_PROD_EV) 
491: ! Fix number of zeros 
492:          NUM_ZERO_EVS = 6 
493: ! Symmetry and rotational component 
494:          ITDET=0.0D0 
495:          CALL DETSYMMETRY(P, HORDER, IT, .FALSE.) 
496:          IF (USEROT) CALL MYINERTIA(P,ITDET) 
497:          IF (ITDET .NE. 0.0D0) THEN 
498:             ITDET=LOG(ITDET)/2.0D0 
499:          END IF 
500: ! Printing 
501:          FEBH_POT_ENE=EREAL 
502:          IF (FETEMP .GT. 0.D0) THEN 
503:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'FE temp=', FETEMP 
504:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Potential energy=', FEBH_POT_ENE 
505:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Log prod=', LOG_PROD_EV 
506:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Horder=', HORDER 
507:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Symfctr=', SYMFCTR 
508:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'ITDET=', ITDET 
509:             IF (DEBUG) WRITE(MYUNIT, '(A,I20)') 'NATOMS=', NATOMS 
510:             IF (DEBUG) WRITE(MYUNIT, '(A,I20)') 'NUM_ZERO_EVS=', NUM_ZERO_EVS 
511:             EREAL=-FETEMP*(SYMFCTR+LOG(1.0D0/HORDER) - EREAL/FETEMP + (3*NATOMS-NUM_ZERO_EVS)*DLOG(FETEMP) - LOG_PROD_EV + ITDET) 
512:             IF (DEBUG) WRITE(MYUNIT, '(A,F20.12)') 'Harmonic superposition contribution=', EREAL - FEBH_POT_ENE 
513:          END IF 
514:       ENDIF488:       ENDIF
515: 489: 
516: ! ds656> If (relative) chemical potentials are present, then we tack490: ! ds656> If (relative) chemical potentials are present, then we tack
517: !        them on to the (free) energy for semi-grand canonical (FE)BH.491: !        them on to the (free) energy for semi-grand canonical (FE)BH.
518:       IF (SEMIGRAND_MUT) THEN492:       IF (SEMIGRAND_MUT) THEN
519:          DUMMY=0.0D0493:          DUMMY=0.0D0
520:          DO J1=2, NSPECIES(0)494:          DO J1=2, NSPECIES(0)
521:             DUMMY = DUMMY + NSPECIES(J1)*SEMIGRAND_MU(J1)495:             DUMMY = DUMMY + NSPECIES(J1)*SEMIGRAND_MU(J1)
522:          ENDDO496:          ENDDO
523:          EREAL = EREAL - DUMMY497:          EREAL = EREAL - DUMMY


r31980/shifth.f90 2017-02-27 18:30:21.232111218 +0000 r31979/shifth.f90 2017-02-27 18:30:23.744144995 +0000
 19:     INTEGER                                 :: N, I 19:     INTEGER                                 :: N, I
 20:     INTEGER                                 :: X, Y, Z 20:     INTEGER                                 :: X, Y, Z
 21:     REAL(REAL64), DIMENSION(:), ALLOCATABLE :: MASSES 21:     REAL(REAL64), DIMENSION(:), ALLOCATABLE :: MASSES
 22:     REAL(REAL64), DIMENSION(:), ALLOCATABLE :: COORDS 22:     REAL(REAL64), DIMENSION(:), ALLOCATABLE :: COORDS
 23:     REAL(REAL64), DIMENSION(3)              :: COM 23:     REAL(REAL64), DIMENSION(3)              :: COM
 24:     REAL(REAL64), DIMENSION(3, 3)           :: INERTIA_TENS 24:     REAL(REAL64), DIMENSION(3, 3)           :: INERTIA_TENS
 25:     REAL(REAL64), DIMENSION(3)              :: I_123 25:     REAL(REAL64), DIMENSION(3)              :: I_123
 26:     REAL(REAL64), DIMENSION(3, 3)           :: PRINCIPAL_AXES 26:     REAL(REAL64), DIMENSION(3, 3)           :: PRINCIPAL_AXES
 27:     REAL(REAL64)                            :: M 27:     REAL(REAL64)                            :: M
 28: ! DSYEV variables 28: ! DSYEV variables
 29:     REAL(REAL64), DIMENSION(100000)           :: WORK 29:     REAL(REAL64), DIMENSION(1000)           :: WORK
 30:     INTEGER, PARAMETER                      :: LWORK = 100000 30:     INTEGER, PARAMETER                      :: LWORK = 1000
 31:     INTEGER                                 :: INFO 31:     INTEGER                                 :: INFO
 32: ! Parameters 32: ! Parameters
 33:     REAL(REAL64), PARAMETER                 :: MAX_INERTIA_ZEROS = 1.0D-3 33:     REAL(REAL64), PARAMETER                 :: MAX_INERTIA_ZEROS = 1.0D-3
 34:  34: 
 35: ! Check that the number of coordinates and size of hessian correspond 35: ! Check that the number of coordinates and size of hessian correspond
 36:     NUM_COORDS = SIZE(INPUT_COORDS) 36:     NUM_COORDS = SIZE(INPUT_COORDS)
 37:     IF (NUM_COORDS .NE. SIZE(HESS, 1)) THEN 37:     IF (NUM_COORDS .NE. SIZE(HESS, 1)) THEN
 38:         STOP 'Coords/hessian size mismatch in subroutine SHIFT_HESS_ZEROS.' 38:         STOP 'Coords/hessian size mismatch in subroutine SHIFT_HESS_ZEROS.'
 39:     END IF 39:     END IF
 40:  40: 


r31980/sparse_bench.F90 2017-02-27 18:30:21.452114176 +0000 r31979/sparse_bench.F90 2017-02-27 18:30:23.968148007 +0000
  1: SUBROUTINE BENCH_SPARSE(COORDINATES, QUENCH_NUM, HESS_CUTOFF, POT_ENERGY, LOG_PROD)  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/sparse_bench.F90' in revision 31979
  2:    USE COMMONS, ONLY: NATOMS, NOFILTER, DEBUG 
  3:    USE MODHESS 
  4: #ifdef __SPARSE 
  5:    USE MODSPARSEHESS 
  6:    USE SHIFT_HESS 
  7: #endif 
  8:    IMPLICIT NONE 
  9: ! Arguments 
 10:    DOUBLE PRECISION, INTENT(IN)  :: COORDINATES(3*NATOMS) 
 11:    INTEGER, INTENT(IN)           :: QUENCH_NUM 
 12:    DOUBLE PRECISION, INTENT(IN)  :: HESS_CUTOFF 
 13:    DOUBLE PRECISION, INTENT(OUT) :: POT_ENERGY 
 14:    DOUBLE PRECISION, INTENT(OUT) :: LOG_PROD 
 15: ! Timing variables 
 16:    DOUBLE PRECISION              :: TIME_DSYEV, TIME_FILT, TIME_SPARSE 
 17:    DOUBLE PRECISION              :: T_START, T_END 
 18: ! File variables 
 19:    INTEGER                       :: SUMMARY_FILE, EFILE1, EFILE2, EFILE3 
 20:    INTEGER                       :: GETUNIT 
 21: ! DSYEV variables 
 22:    INTEGER, PARAMETER            :: LWORK = 100000 
 23:    DOUBLE PRECISION              :: WORK(LWORK) 
 24:    INTEGER                       :: INFO 
 25: ! Other variables 
 26:    DOUBLE PRECISION              :: LOG_PROD_DSYEV, LOG_PROD_FILT, LOG_PROD_SPARSE 
 27:    DOUBLE PRECISION              :: EVALUES(3*NATOMS) 
 28:    DOUBLE PRECISION              :: SHIFT_ARRAY(6) 
 29:    DOUBLE PRECISION              :: GRAD(3*NATOMS) 
 30:    DOUBLE PRECISION, ALLOCATABLE :: HESS_COPY(:, :), COORDSCOPY(:) 
 31:    INTEGER                       :: NUM_ZERO_EVS , MYVAR, I, J 
 32:  
 33:    INTEGER, PARAMETER              :: HESS_OUT1 = 4739 
 34:  
 35: #ifdef __SPARSE 
 36: ! Assume 6 zero eigenvalues 
 37:    NUM_ZERO_EVS = 6 
 38:  
 39: ! Calculate hessian and mass weight 
 40:    IF (ALLOCATED(HESS)) DEALLOCATE(HESS) 
 41:    ALLOCATE(HESS(3*NATOMS, 3*NATOMS)) 
 42:    CALL POTENTIAL(COORDINATES, GRAD, POT_ENERGY, .TRUE., .TRUE.) 
 43:    CALL MASSWT() 
 44:  
 45: ! Copy the original (unfiltered) hessian to somewhere safe 
 46:    IF (ALLOCATED(HESS_COPY)) DEALLOCATE(HESS_COPY) 
 47:    ALLOCATE(HESS_COPY(3*NATOMS, 3*NATOMS)) 
 48:    HESS_COPY = HESS 
 49:  
 50: ! Copy the original coordinates to somewhere safe 
 51:    IF (ALLOCATED(COORDSCOPY)) DEALLOCATE(COORDSCOPY) 
 52:    ALLOCATE(COORDSCOPY(3*NATOMS)) 
 53:    COORDSCOPY = COORDINATES 
 54:  
 55: ! Assign shift array 
 56:    SHIFT_ARRAY(1:6) = 1.0D0 
 57:  
 58: ! Sparse section 
 59:    CALL CPU_TIME(T_START) 
 60:    CALL SHIFT_HESS_ZEROS(COORDINATES, SHIFT_ARRAY) 
 61:    CALL FILTER_ZEROS(HESS, HESS_CUTOFF) 
 62:    CALL GET_DETERMINANT(3*NATOMS, LOG_PROD_SPARSE, QUENCH_NUM) 
 63:    CALL CPU_TIME(T_END) 
 64:    TIME_SPARSE = T_END - T_START 
 65:  
 66: ! DSYEV without filters section 
 67:    HESS = HESS_COPY 
 68:    CALL CPU_TIME(T_START) 
 69:    CALL SHIFT_HESS_ZEROS(COORDINATES, SHIFT_ARRAY) 
 70:    CALL DSYEV('N', 'L', 3*NATOMS, HESS, 3*NATOMS, EVALUES, WORK, LWORK, INFO) 
 71:    CALL CPU_TIME(T_END) 
 72:    TIME_DSYEV = T_END - T_START 
 73:    LOG_PROD_DSYEV = SUM(DLOG(EVALUES(NUM_ZERO_EVS + 1:))) 
 74:  
 75: ! Need to pass back the log product of frequencies for actual FEBH (sqrt of eigenvalues) 
 76:    LOG_PROD = 0.5D0 * LOG_PROD_DSYEV 
 77:  
 78: ! Print summary 
 79:    SUMMARY_FILE = GETUNIT()  
 80:    OPEN(UNIT=SUMMARY_FILE, FILE='sparse_benchmarks', STATUS='UNKNOWN', ACCESS='APPEND') 
 81:    WRITE(SUMMARY_FILE, '(A, I12)') 'Quench', QUENCH_NUM 
 82:    WRITE(SUMMARY_FILE, '(A, A, E18.10, A, F10.4, A, E18.10, A, F8.2)') 'DSYEV             !!||', & 
 83:                                                                     ' log_prod= ', LOG_PROD_DSYEV, & 
 84:                                                                     ' time= ', TIME_DSYEV, & 
 85:                                                                     ' error= ', LOG_PROD_DSYEV - LOG_PROD_DSYEV, & 
 86:                                                                     ' % error= ', 100.0D0 * (LOG_PROD_DSYEV - LOG_PROD_DSYEV) / LOG_PROD_DSYEV 
 87:    WRITE(SUMMARY_FILE, '(A, A, E18.10, A, F10.4, A, E18.10, A, F8.2)') 'Sparse            ||', & 
 88:                                                                     ' log_prod= ', LOG_PROD_SPARSE, & 
 89:                                                                     ' time= ', TIME_SPARSE, & 
 90:                                                                     ' error= ', LOG_PROD_SPARSE - LOG_PROD_DSYEV, & 
 91:                                                                     ' % error= ', 100.0D0 * (LOG_PROD_SPARSE - LOG_PROD_DSYEV) / LOG_PROD_DSYEV 
 92:    CLOSE(SUMMARY_FILE)  
 93: #endif  
 94:  
 95: END SUBROUTINE BENCH_SPARSE 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0