add_subdirectory(source/lib/xdrf)
-if(UNRES_NA_MMCE)
+if(UNRES_NA_MMCE) #kompiluj na mmce
if(UNRES_WITH_MPI)
# Brak MPI dla gfortrana, wiec tylko na ifort sie skompiluje
add_subdirectory(source/xdrfpdb/src)
add_subdirectory(source/xdrfpdb/src-M)
add_subdirectory(source/maxlik/src_CSA)
-else()
+else() #kompiluj gdzie indziej
add_subdirectory(source/unres/src_MD)
if(UNRES_WITH_MPI)
parmread.F
pinorm.f
printmat.f
- prng_32.F
+ prng.F
randgens.f
ran.f
readpdb.F
-Makefile-DFA-NEWPARM.piasek
\ No newline at end of file
+Makefile_4P
\ No newline at end of file
CPPFLAGS = -DPROCOR -DLINUX -DPGI -DISNAN -DMP -DMPI -DUNRES \
-DSPLITELE -DAMD64 -DLANG0 -DPROCOR \
- -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
-# -DPROCOR
+ -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DPROCOR
# -DTSCSC
#-DTIMING \
# -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
distfit.o banach.o TMscore_subroutine.o minim_mult.o
unres: ${object}
-# cc -o compinfo compinfo.c
-# ./compinfo | true
+ cc -o compinfo compinfo.c
+ ./compinfo | true
${FC} ${FFLAGS} cinfo.f
${FC} ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN}
--- /dev/null
+C 0 4196011 1
+ subroutine cinfo
+ include 'COMMON.IOUNITS'
+ write(iout,*)'++++ Compile info ++++'
+ write(iout,*)'Version 0.4196011 build 1'
+ write(iout,*)'compiled Mon Nov 5 08:27:57 2012'
+ write(iout,*)'compiled by bz94@matrix.chem.cornell.edu'
+ write(iout,*)'OS name: Linux '
+ write(iout,*)'OS release: 2.6.34.9-69.fc13.x86_64 '
+ write(iout,*)'OS version:',
+ & ' #1 SMP Tue May 3 09:23:03 UTC 2011 '
+ write(iout,*)'flags:'
+ write(iout,*)'CPPFLAGS = -DPROCOR -DLINUX -DPGI -DISNAN -DMP ...'
+ write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...'
+ write(iout,*)'FC= ifort '
+ write(iout,*)'OPT = -O3 -ip -w '
+ write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include'
+ write(iout,*)'FFLAGS1 = -c -w -g -O0 -d2 -CA -CB -I$(INSTALL_...'
+ write(iout,*)'FFLAGS2 = -c -w -g -O0 -I$(INSTALL_DIR)/include'
+ write(iout,*)'FFLAGS3 = -c -w -O3 -mp'
+ write(iout,*)'FFLAGSE = -c -w -O3 -ipo -ipo_obj -opt_report ...'
+ write(iout,*)'BIN = ../../../bin/unres/CSA/unres_csa_ifort_mp...'
+ write(iout,*)'LIBS = -lpthread -L$(INSTALL_DIR)/lib -lmpich'
+ write(iout,*)'ARCH = LINUX'
+ write(iout,*)'PP = /lib/cpp -P'
+ write(iout,*)'object = unres_csa.o arcos.o cartprint.o chainb...'
+ write(iout,*)'++++ End of compile info ++++'
+ return
+ end
--- /dev/null
+ real*8 function prng_next(me)
+ implicit none
+ integer me
+c
+c Calling sequence:
+c <new random number> = prng_next ( <ordinal of generator desired> )
+c <vector of random #s> = vprng ( <ordinal>, <vector>, <length> )
+c
+c This code is based on a sequential algorithm provided by Mal Kalos.
+c This version uses a single 64-bit word to store the initial seeds
+c and additive constants.
+c A 64-bit floating point number is returned.
+c
+c The array "iparam" is full-word aligned, being padded by zeros to
+c let each generator be on a subpage boundary.
+c That is, rows 1 and 2 in a given column of the array are for real,
+c rows 3-16 are bogus.
+c
+c July 12, 1993: double the number of sequences. We should have been
+c using two packets per seed, rather than four
+c October 31, 1993: merge the two arrays of seeds and constants,
+c and switch to 64-bit arithmetic.
+c June 1994: port to RS6K. Internal state is kept as 2 64-bit integers
+c The ishft function is defined only on 32-bit integers, so we will
+c shift numbers by dividing by 2**11 and then adding on 2**53-1.
+c
+c November 1994: ishift now works on 64-bit numbers (though it gives a
+c warning). Thus we go back to using it. John Zollweg also added the
+c vprng() routine to return vectors of real*8 random numbers.
+c
+ real*8 recip53
+ parameter ( recip53 = 2.0D0**(-53) )
+ integer*8 two
+ parameter ( two = 2**11)
+ integer*8 m,ishift
+c parameter ( m = 34522712143931 ) ! 11**13
+c parameter ( ishift = 9007199254740991 ) ! 2**53-1
+
+ integer nmax
+ integer*8 iparam
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ integer*8 next
+
+crc g77 doesn't support integer*8 constants
+ m = dint(34522712143931.0d0)
+ ishift = dint(9007199254740991.0d0)
+
+c RS6K porting note: ishift now takes 64-bit integers , with a warning
+ if ( 0.le.me .and. me.le.nmax ) then
+ next = iparam(1,me)*m + iparam(2,me)
+ iparam(1,me) = next
+ prng_next = recip53 * ishft( next, -11 )
+ else
+ prng_next=-1.0D0
+ endif
+
+ end
+c
+c vprng(me, rn, num) Get a vector of random numbers
+c
+ subroutine vprng(me,rn,num)
+ real*8 recip53, rn(1)
+ parameter ( recip53 = 2.0D0**(-53) )
+ integer*8 m,iparam
+c parameter ( m = 34522712143931 ) ! 11**13
+ integer nmax, num, me
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ integer*8 next
+
+crc g77 doesn't support integer*8 constants
+ m = dint(34522712143931.0d0)
+
+ if ( 0.le.me .and. me.le.nmax ) then
+ do 1 i=1,num
+ next = iparam(1,me)*m + iparam(2,me)
+ iparam(1,me) = next
+ rn(i) = recip53 * ishft( next, -11 )
+ 1 continue
+ else
+ rn(1)=-1.0D0
+ endif
+ return
+ end
+
+c
+c prng_chkpnt Get the current state of a generator
+c
+c Calling sequence:
+c logical prng_chkpnt, status
+c status = prng_chkpnt (me, iseed) where
+c
+c me is the particular generator whose state is being gotten
+c seed is an 4-element integer array where the "l"-values will be saved
+c
+ logical function prng_chkpnt (me, iseed)
+ implicit none
+ integer me
+ integer*8 iseed
+
+ integer nmax
+ integer*8 iparam
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ if (me .lt. 0 .or. me .gt. nmax) then
+ prng_chkpnt=.false.
+ else
+ prng_chkpnt=.true.
+ iseed=iparam(1,me)
+ endif
+ end
+c
+c prng_restart Restart generator from a saved state
+c
+c Calling sequence:
+c logical prng_restart, status
+c status = prng_restart (me, iseed) where
+c
+c me is the particular generator being restarted
+c iseed is a 8-byte integer containing the "l"-values
+c
+ logical function prng_restart (me, iseed)
+ implicit none
+ integer me
+ integer*8 iseed
+
+ integer nmax
+ integer*8 iparam
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ if (me .lt. 0 .or. me .gt. nmax) then
+ prng_restart=.false.
+ return
+ else
+ prng_restart=.true.
+ iparam(1,me)=iseed
+ endif
+ end
+
+ block data prngblk
+ parameter(nmax=1021)
+ integer*8 iparam
+ common/ksrprng/iparam(2,0:nmax)
+ data (iparam(1,i),iparam(2,i),i= 0, 29) /
+ + 11848219, 11848219, 11848237, 11848237, 11848241, 11848241,
+ + 11848247, 11848247, 11848253, 11848253, 11848271, 11848271,
+ + 11848297, 11848297, 11848313, 11848313, 11848339, 11848339,
+ + 11848351, 11848351, 11848357, 11848357, 11848363, 11848363,
+ + 11848367, 11848367, 11848373, 11848373, 11848379, 11848379,
+ + 11848393, 11848393, 11848433, 11848433, 11848451, 11848451,
+ + 11848469, 11848469, 11848477, 11848477, 11848489, 11848489,
+ + 11848493, 11848493, 11848513, 11848513, 11848523, 11848523,
+ + 11848531, 11848531, 11848537, 11848537, 11848553, 11848553,
+ + 11848589, 11848589, 11848591, 11848591, 11848601, 11848601 /
+ data (iparam(1,i),iparam(2,i),i= 30, 59) /
+ + 11848619, 11848619, 11848637, 11848637, 11848663, 11848663,
+ + 11848673, 11848673, 11848679, 11848679, 11848691, 11848691,
+ + 11848699, 11848699, 11848709, 11848709, 11848717, 11848717,
+ + 11848721, 11848721, 11848729, 11848729, 11848741, 11848741,
+ + 11848751, 11848751, 11848757, 11848757, 11848787, 11848787,
+ + 11848801, 11848801, 11848829, 11848829, 11848853, 11848853,
+ + 11848861, 11848861, 11848867, 11848867, 11848873, 11848873,
+ + 11848891, 11848891, 11848909, 11848909, 11848919, 11848919,
+ + 11848931, 11848931, 11848937, 11848937, 11848961, 11848961,
+ + 11848981, 11848981, 11849021, 11849021, 11849039, 11849039 /
+ data (iparam(1,i),iparam(2,i),i= 60, 89) /
+ + 11849053, 11849053, 11849059, 11849059, 11849069, 11849069,
+ + 11849077, 11849077, 11849087, 11849087, 11849093, 11849093,
+ + 11849107, 11849107, 11849111, 11849111, 11849129, 11849129,
+ + 11849137, 11849137, 11849177, 11849177, 11849183, 11849183,
+ + 11849203, 11849203, 11849231, 11849231, 11849237, 11849237,
+ + 11849239, 11849239, 11849249, 11849249, 11849251, 11849251,
+ + 11849269, 11849269, 11849273, 11849273, 11849291, 11849291,
+ + 11849297, 11849297, 11849309, 11849309, 11849339, 11849339,
+ + 11849359, 11849359, 11849363, 11849363, 11849399, 11849399,
+ + 11849401, 11849401, 11849413, 11849413, 11849417, 11849417 /
+ data (iparam(1,i),iparam(2,i),i= 90, 119) /
+ + 11849437, 11849437, 11849443, 11849443, 11849473, 11849473,
+ + 11849491, 11849491, 11849503, 11849503, 11849507, 11849507,
+ + 11849557, 11849557, 11849567, 11849567, 11849569, 11849569,
+ + 11849573, 11849573, 11849587, 11849587, 11849599, 11849599,
+ + 11849633, 11849633, 11849641, 11849641, 11849653, 11849653,
+ + 11849659, 11849659, 11849671, 11849671, 11849683, 11849683,
+ + 11849689, 11849689, 11849693, 11849693, 11849699, 11849699,
+ + 11849701, 11849701, 11849707, 11849707, 11849713, 11849713,
+ + 11849723, 11849723, 11849741, 11849741, 11849743, 11849743,
+ + 11849759, 11849759, 11849767, 11849767, 11849771, 11849771 /
+ data (iparam(1,i),iparam(2,i),i= 120, 149) /
+ + 11849791, 11849791, 11849801, 11849801, 11849809, 11849809,
+ + 11849813, 11849813, 11849869, 11849869, 11849881, 11849881,
+ + 11849891, 11849891, 11849909, 11849909, 11849923, 11849923,
+ + 11849933, 11849933, 11849947, 11849947, 11849987, 11849987,
+ + 11850001, 11850001, 11850011, 11850011, 11850019, 11850019,
+ + 11850023, 11850023, 11850031, 11850031, 11850049, 11850049,
+ + 11850061, 11850061, 11850073, 11850073, 11850077, 11850077,
+ + 11850103, 11850103, 11850109, 11850109, 11850121, 11850121,
+ + 11850127, 11850127, 11850133, 11850133, 11850149, 11850149,
+ + 11850161, 11850161, 11850169, 11850169, 11850191, 11850191 /
+ data (iparam(1,i),iparam(2,i),i= 150, 179) /
+ + 11850233, 11850233, 11850247, 11850247, 11850259, 11850259,
+ + 11850269, 11850269, 11850283, 11850283, 11850301, 11850301,
+ + 11850341, 11850341, 11850347, 11850347, 11850367, 11850367,
+ + 11850373, 11850373, 11850379, 11850379, 11850389, 11850389,
+ + 11850407, 11850407, 11850427, 11850427, 11850437, 11850437,
+ + 11850469, 11850469, 11850481, 11850481, 11850511, 11850511,
+ + 11850529, 11850529, 11850541, 11850541, 11850557, 11850557,
+ + 11850607, 11850607, 11850611, 11850611, 11850667, 11850667,
+ + 11850677, 11850677, 11850679, 11850679, 11850701, 11850701,
+ + 11850731, 11850731, 11850739, 11850739, 11850749, 11850749 /
+ data (iparam(1,i),iparam(2,i),i= 180, 209) /
+ + 11850791, 11850791, 11850803, 11850803, 11850829, 11850829,
+ + 11850833, 11850833, 11850859, 11850859, 11850877, 11850877,
+ + 11850899, 11850899, 11850907, 11850907, 11850913, 11850913,
+ + 11850919, 11850919, 11850931, 11850931, 11850941, 11850941,
+ + 11850947, 11850947, 11850953, 11850953, 11850961, 11850961,
+ + 11850983, 11850983, 11850991, 11850991, 11850997, 11850997,
+ + 11851031, 11851031, 11851033, 11851033, 11851051, 11851051,
+ + 11851061, 11851061, 11851067, 11851067, 11851093, 11851093,
+ + 11851109, 11851109, 11851123, 11851123, 11851127, 11851127,
+ + 11851139, 11851139, 11851157, 11851157, 11851163, 11851163 /
+ data (iparam(1,i),iparam(2,i),i= 210, 239) /
+ + 11851181, 11851181, 11851201, 11851201, 11851219, 11851219,
+ + 11851291, 11851291, 11851303, 11851303, 11851309, 11851309,
+ + 11851313, 11851313, 11851319, 11851319, 11851349, 11851349,
+ + 11851351, 11851351, 11851361, 11851361, 11851373, 11851373,
+ + 11851403, 11851403, 11851409, 11851409, 11851423, 11851423,
+ + 11851447, 11851447, 11851451, 11851451, 11851481, 11851481,
+ + 11851493, 11851493, 11851519, 11851519, 11851523, 11851523,
+ + 11851529, 11851529, 11851547, 11851547, 11851549, 11851549,
+ + 11851559, 11851559, 11851577, 11851577, 11851589, 11851589,
+ + 11851591, 11851591, 11851597, 11851597, 11851603, 11851603 /
+ data (iparam(1,i),iparam(2,i),i= 240, 269) /
+ + 11851607, 11851607, 11851613, 11851613, 11851621, 11851621,
+ + 11851627, 11851627, 11851639, 11851639, 11851673, 11851673,
+ + 11851681, 11851681, 11851727, 11851727, 11851753, 11851753,
+ + 11851759, 11851759, 11851787, 11851787, 11851793, 11851793,
+ + 11851799, 11851799, 11851813, 11851813, 11851841, 11851841,
+ + 11851859, 11851859, 11851867, 11851867, 11851891, 11851891,
+ + 11851909, 11851909, 11851919, 11851919, 11851927, 11851927,
+ + 11851933, 11851933, 11851949, 11851949, 11851967, 11851967,
+ + 11851997, 11851997, 11852017, 11852017, 11852051, 11852051,
+ + 11852053, 11852053, 11852059, 11852059, 11852083, 11852083 /
+ data (iparam(1,i),iparam(2,i),i= 270, 299) /
+ + 11852089, 11852089, 11852129, 11852129, 11852147, 11852147,
+ + 11852149, 11852149, 11852161, 11852161, 11852171, 11852171,
+ + 11852177, 11852177, 11852209, 11852209, 11852221, 11852221,
+ + 11852237, 11852237, 11852251, 11852251, 11852263, 11852263,
+ + 11852273, 11852273, 11852279, 11852279, 11852287, 11852287,
+ + 11852293, 11852293, 11852297, 11852297, 11852303, 11852303,
+ + 11852311, 11852311, 11852327, 11852327, 11852339, 11852339,
+ + 11852341, 11852341, 11852359, 11852359, 11852369, 11852369,
+ + 11852437, 11852437, 11852453, 11852453, 11852459, 11852459,
+ + 11852473, 11852473, 11852513, 11852513, 11852531, 11852531 /
+ data (iparam(1,i),iparam(2,i),i= 300, 329) /
+ + 11852537, 11852537, 11852539, 11852539, 11852557, 11852557,
+ + 11852573, 11852573, 11852579, 11852579, 11852591, 11852591,
+ + 11852609, 11852609, 11852611, 11852611, 11852623, 11852623,
+ + 11852641, 11852641, 11852647, 11852647, 11852657, 11852657,
+ + 11852663, 11852663, 11852717, 11852717, 11852719, 11852719,
+ + 11852741, 11852741, 11852759, 11852759, 11852767, 11852767,
+ + 11852773, 11852773, 11852803, 11852803, 11852807, 11852807,
+ + 11852809, 11852809, 11852831, 11852831, 11852833, 11852833,
+ + 11852837, 11852837, 11852857, 11852857, 11852873, 11852873,
+ + 11852879, 11852879, 11852891, 11852891, 11852917, 11852917 /
+ data (iparam(1,i),iparam(2,i),i= 330, 359) /
+ + 11852921, 11852921, 11852957, 11852957, 11852959, 11852959,
+ + 11852969, 11852969, 11852983, 11852983, 11852989, 11852989,
+ + 11853001, 11853001, 11853013, 11853013, 11853019, 11853019,
+ + 11853031, 11853031, 11853089, 11853089, 11853133, 11853133,
+ + 11853157, 11853157, 11853161, 11853161, 11853181, 11853181,
+ + 11853203, 11853203, 11853217, 11853217, 11853221, 11853221,
+ + 11853227, 11853227, 11853241, 11853241, 11853307, 11853307,
+ + 11853319, 11853319, 11853323, 11853323, 11853329, 11853329,
+ + 11853367, 11853367, 11853383, 11853383, 11853419, 11853419,
+ + 11853421, 11853421, 11853427, 11853427, 11853449, 11853449 /
+ data (iparam(1,i),iparam(2,i),i= 360, 389) /
+ + 11853451, 11853451, 11853463, 11853463, 11853529, 11853529,
+ + 11853557, 11853557, 11853571, 11853571, 11853601, 11853601,
+ + 11853613, 11853613, 11853617, 11853617, 11853629, 11853629,
+ + 11853649, 11853649, 11853659, 11853659, 11853679, 11853679,
+ + 11853689, 11853689, 11853719, 11853719, 11853731, 11853731,
+ + 11853757, 11853757, 11853761, 11853761, 11853773, 11853773,
+ + 11853791, 11853791, 11853817, 11853817, 11853839, 11853839,
+ + 11853847, 11853847, 11853857, 11853857, 11853869, 11853869,
+ + 11853883, 11853883, 11853887, 11853887, 11853889, 11853889,
+ + 11853893, 11853893, 11853899, 11853899, 11853911, 11853911 /
+ data (iparam(1,i),iparam(2,i),i= 390, 419) /
+ + 11853931, 11853931, 11853943, 11853943, 11853979, 11853979,
+ + 11853991, 11853991, 11854001, 11854001, 11854009, 11854009,
+ + 11854019, 11854019, 11854057, 11854057, 11854061, 11854061,
+ + 11854147, 11854147, 11854159, 11854159, 11854163, 11854163,
+ + 11854169, 11854169, 11854211, 11854211, 11854247, 11854247,
+ + 11854261, 11854261, 11854267, 11854267, 11854279, 11854279,
+ + 11854303, 11854303, 11854327, 11854327, 11854331, 11854331,
+ + 11854333, 11854333, 11854363, 11854363, 11854379, 11854379,
+ + 11854399, 11854399, 11854411, 11854411, 11854429, 11854429,
+ + 11854433, 11854433, 11854439, 11854439, 11854441, 11854441 /
+ data (iparam(1,i),iparam(2,i),i= 420, 449) /
+ + 11854463, 11854463, 11854477, 11854477, 11854489, 11854489,
+ + 11854517, 11854517, 11854519, 11854519, 11854523, 11854523,
+ + 11854529, 11854529, 11854567, 11854567, 11854571, 11854571,
+ + 11854573, 11854573, 11854603, 11854603, 11854607, 11854607,
+ + 11854681, 11854681, 11854691, 11854691, 11854709, 11854709,
+ + 11854723, 11854723, 11854757, 11854757, 11854783, 11854783,
+ + 11854793, 11854793, 11854813, 11854813, 11854847, 11854847,
+ + 11854853, 11854853, 11854873, 11854873, 11854877, 11854877,
+ + 11854883, 11854883, 11854891, 11854891, 11854897, 11854897,
+ + 11854901, 11854901, 11854919, 11854919, 11854937, 11854937 /
+ data (iparam(1,i),iparam(2,i),i= 450, 479) /
+ + 11854961, 11854961, 11854963, 11854963, 11854979, 11854979,
+ + 11855003, 11855003, 11855017, 11855017, 11855023, 11855023,
+ + 11855029, 11855029, 11855033, 11855033, 11855111, 11855111,
+ + 11855141, 11855141, 11855147, 11855147, 11855149, 11855149,
+ + 11855159, 11855159, 11855177, 11855177, 11855203, 11855203,
+ + 11855213, 11855213, 11855219, 11855219, 11855231, 11855231,
+ + 11855267, 11855267, 11855269, 11855269, 11855303, 11855303,
+ + 11855309, 11855309, 11855321, 11855321, 11855329, 11855329,
+ + 11855339, 11855339, 11855351, 11855351, 11855353, 11855353,
+ + 11855357, 11855357, 11855359, 11855359, 11855381, 11855381 /
+ data (iparam(1,i),iparam(2,i),i= 480, 509) /
+ + 11855383, 11855383, 11855387, 11855387, 11855399, 11855399,
+ + 11855407, 11855407, 11855413, 11855413, 11855489, 11855489,
+ + 11855491, 11855491, 11855507, 11855507, 11855521, 11855521,
+ + 11855531, 11855531, 11855549, 11855549, 11855551, 11855551,
+ + 11855567, 11855567, 11855581, 11855581, 11855587, 11855587,
+ + 11855593, 11855593, 11855633, 11855633, 11855653, 11855653,
+ + 11855663, 11855663, 11855687, 11855687, 11855689, 11855689,
+ + 11855699, 11855699, 11855713, 11855713, 11855731, 11855731,
+ + 11855737, 11855737, 11855743, 11855743, 11855747, 11855747,
+ + 11855759, 11855759, 11855773, 11855773, 11855801, 11855801 /
+ data (iparam(1,i),iparam(2,i),i= 510, 539) /
+ + 11855807, 11855807, 11855813, 11855813, 11855827, 11855827,
+ + 11855839, 11855839, 11855869, 11855869, 11855881, 11855881,
+ + 11855903, 11855903, 11855911, 11855911, 11855933, 11855933,
+ + 11855959, 11855959, 11855989, 11855989, 11855993, 11855993,
+ + 11855999, 11855999, 11856001, 11856001, 11856023, 11856023,
+ + 11856049, 11856049, 11856071, 11856071, 11856101, 11856101,
+ + 11856107, 11856107, 11856113, 11856113, 11856139, 11856139,
+ + 11856151, 11856151, 11856161, 11856161, 11856179, 11856179,
+ + 11856193, 11856193, 11856199, 11856199, 11856223, 11856223,
+ + 11856239, 11856239, 11856263, 11856263, 11856269, 11856269 /
+ data (iparam(1,i),iparam(2,i),i= 540, 569) /
+ + 11856281, 11856281, 11856287, 11856287, 11856307, 11856307,
+ + 11856311, 11856311, 11856329, 11856329, 11856343, 11856343,
+ + 11856359, 11856359, 11856371, 11856371, 11856373, 11856373,
+ + 11856409, 11856409, 11856419, 11856419, 11856461, 11856461,
+ + 11856469, 11856469, 11856473, 11856473, 11856479, 11856479,
+ + 11856511, 11856511, 11856517, 11856517, 11856541, 11856541,
+ + 11856547, 11856547, 11856553, 11856553, 11856583, 11856583,
+ + 11856629, 11856629, 11856641, 11856641, 11856653, 11856653,
+ + 11856659, 11856659, 11856673, 11856673, 11856697, 11856697,
+ + 11856709, 11856709, 11856727, 11856727, 11856731, 11856731 /
+ data (iparam(1,i),iparam(2,i),i= 570, 599) /
+ + 11856763, 11856763, 11856809, 11856809, 11856811, 11856811,
+ + 11856821, 11856821, 11856841, 11856841, 11856857, 11856857,
+ + 11856877, 11856877, 11856883, 11856883, 11856899, 11856899,
+ + 11856919, 11856919, 11856947, 11856947, 11856953, 11856953,
+ + 11856979, 11856979, 11857003, 11857003, 11857033, 11857033,
+ + 11857037, 11857037, 11857039, 11857039, 11857049, 11857049,
+ + 11857061, 11857061, 11857067, 11857067, 11857073, 11857073,
+ + 11857081, 11857081, 11857091, 11857091, 11857093, 11857093,
+ + 11857099, 11857099, 11857123, 11857123, 11857127, 11857127,
+ + 11857147, 11857147, 11857151, 11857151, 11857193, 11857193 /
+ data (iparam(1,i),iparam(2,i),i= 600, 629) /
+ + 11857217, 11857217, 11857229, 11857229, 11857243, 11857243,
+ + 11857249, 11857249, 11857267, 11857267, 11857277, 11857277,
+ + 11857291, 11857291, 11857303, 11857303, 11857309, 11857309,
+ + 11857327, 11857327, 11857331, 11857331, 11857333, 11857333,
+ + 11857361, 11857361, 11857367, 11857367, 11857369, 11857369,
+ + 11857393, 11857393, 11857399, 11857399, 11857409, 11857409,
+ + 11857421, 11857421, 11857423, 11857423, 11857451, 11857451,
+ + 11857453, 11857453, 11857457, 11857457, 11857477, 11857477,
+ + 11857481, 11857481, 11857493, 11857493, 11857499, 11857499,
+ + 11857519, 11857519, 11857523, 11857523, 11857529, 11857529 /
+ data (iparam(1,i),iparam(2,i),i= 630, 659) /
+ + 11857543, 11857543, 11857561, 11857561, 11857589, 11857589,
+ + 11857591, 11857591, 11857613, 11857613, 11857621, 11857621,
+ + 11857661, 11857661, 11857667, 11857667, 11857693, 11857693,
+ + 11857697, 11857697, 11857709, 11857709, 11857711, 11857711,
+ + 11857751, 11857751, 11857753, 11857753, 11857759, 11857759,
+ + 11857763, 11857763, 11857777, 11857777, 11857787, 11857787,
+ + 11857793, 11857793, 11857801, 11857801, 11857817, 11857817,
+ + 11857819, 11857819, 11857831, 11857831, 11857837, 11857837,
+ + 11857873, 11857873, 11857877, 11857877, 11857883, 11857883,
+ + 11857889, 11857889, 11857907, 11857907, 11857913, 11857913 /
+ data (iparam(1,i),iparam(2,i),i= 660, 689) /
+ + 11857931, 11857931, 11857969, 11857969, 11857991, 11857991,
+ + 11857999, 11857999, 11858009, 11858009, 11858017, 11858017,
+ + 11858023, 11858023, 11858029, 11858029, 11858039, 11858039,
+ + 11858051, 11858051, 11858057, 11858057, 11858059, 11858059,
+ + 11858101, 11858101, 11858111, 11858111, 11858131, 11858131,
+ + 11858149, 11858149, 11858159, 11858159, 11858177, 11858177,
+ + 11858191, 11858191, 11858201, 11858201, 11858227, 11858227,
+ + 11858243, 11858243, 11858267, 11858267, 11858269, 11858269,
+ + 11858279, 11858279, 11858281, 11858281, 11858291, 11858291,
+ + 11858311, 11858311, 11858323, 11858323, 11858359, 11858359 /
+ data (iparam(1,i),iparam(2,i),i= 690, 719) /
+ + 11858377, 11858377, 11858381, 11858381, 11858387, 11858387,
+ + 11858423, 11858423, 11858443, 11858443, 11858447, 11858447,
+ + 11858479, 11858479, 11858533, 11858533, 11858543, 11858543,
+ + 11858551, 11858551, 11858557, 11858557, 11858569, 11858569,
+ + 11858573, 11858573, 11858579, 11858579, 11858597, 11858597,
+ + 11858599, 11858599, 11858629, 11858629, 11858657, 11858657,
+ + 11858659, 11858659, 11858683, 11858683, 11858701, 11858701,
+ + 11858719, 11858719, 11858723, 11858723, 11858729, 11858729,
+ + 11858747, 11858747, 11858779, 11858779, 11858783, 11858783,
+ + 11858801, 11858801, 11858807, 11858807, 11858813, 11858813 /
+ data (iparam(1,i),iparam(2,i),i= 720, 749) /
+ + 11858839, 11858839, 11858851, 11858851, 11858893, 11858893,
+ + 11858897, 11858897, 11858921, 11858921, 11858947, 11858947,
+ + 11858953, 11858953, 11858969, 11858969, 11858971, 11858971,
+ + 11858989, 11858989, 11859017, 11859017, 11859031, 11859031,
+ + 11859049, 11859049, 11859061, 11859061, 11859073, 11859073,
+ + 11859077, 11859077, 11859079, 11859079, 11859083, 11859083,
+ + 11859101, 11859101, 11859109, 11859109, 11859137, 11859137,
+ + 11859139, 11859139, 11859151, 11859151, 11859157, 11859157,
+ + 11859163, 11859163, 11859167, 11859167, 11859179, 11859179,
+ + 11859187, 11859187, 11859229, 11859229, 11859233, 11859233 /
+ data (iparam(1,i),iparam(2,i),i= 750, 779) /
+ + 11859241, 11859241, 11859247, 11859247, 11859269, 11859269,
+ + 11859293, 11859293, 11859307, 11859307, 11859311, 11859311,
+ + 11859349, 11859349, 11859359, 11859359, 11859371, 11859371,
+ + 11859377, 11859377, 11859383, 11859383, 11859427, 11859427,
+ + 11859433, 11859433, 11859451, 11859451, 11859457, 11859457,
+ + 11859461, 11859461, 11859473, 11859473, 11859481, 11859481,
+ + 11859487, 11859487, 11859493, 11859493, 11859503, 11859503,
+ + 11859509, 11859509, 11859539, 11859539, 11859541, 11859541,
+ + 11859563, 11859563, 11859569, 11859569, 11859571, 11859571,
+ + 11859583, 11859583, 11859599, 11859599, 11859611, 11859611 /
+ data (iparam(1,i),iparam(2,i),i= 780, 809) /
+ + 11859643, 11859643, 11859707, 11859707, 11859713, 11859713,
+ + 11859719, 11859719, 11859739, 11859739, 11859751, 11859751,
+ + 11859791, 11859791, 11859817, 11859817, 11859821, 11859821,
+ + 11859833, 11859833, 11859847, 11859847, 11859853, 11859853,
+ + 11859877, 11859877, 11859889, 11859889, 11859893, 11859893,
+ + 11859901, 11859901, 11859907, 11859907, 11859917, 11859917,
+ + 11859923, 11859923, 11859929, 11859929, 11859961, 11859961,
+ + 11859979, 11859979, 11859989, 11859989, 11859997, 11859997,
+ + 11860021, 11860021, 11860031, 11860031, 11860039, 11860039,
+ + 11860049, 11860049, 11860081, 11860081, 11860087, 11860087 /
+ data (iparam(1,i),iparam(2,i),i= 810, 839) /
+ + 11860097, 11860097, 11860103, 11860103, 11860109, 11860109,
+ + 11860117, 11860117, 11860133, 11860133, 11860151, 11860151,
+ + 11860171, 11860171, 11860207, 11860207, 11860223, 11860223,
+ + 11860231, 11860231, 11860243, 11860243, 11860267, 11860267,
+ + 11860301, 11860301, 11860307, 11860307, 11860327, 11860327,
+ + 11860379, 11860379, 11860397, 11860397, 11860411, 11860411,
+ + 11860469, 11860469, 11860477, 11860477, 11860483, 11860483,
+ + 11860487, 11860487, 11860489, 11860489, 11860493, 11860493,
+ + 11860517, 11860517, 11860547, 11860547, 11860567, 11860567,
+ + 11860573, 11860573, 11860613, 11860613, 11860619, 11860619 /
+ data (iparam(1,i),iparam(2,i),i= 840, 869) /
+ + 11860627, 11860627, 11860637, 11860637, 11860643, 11860643,
+ + 11860649, 11860649, 11860661, 11860661, 11860669, 11860669,
+ + 11860687, 11860687, 11860691, 11860691, 11860697, 11860697,
+ + 11860699, 11860699, 11860703, 11860703, 11860727, 11860727,
+ + 11860741, 11860741, 11860753, 11860753, 11860777, 11860777,
+ + 11860787, 11860787, 11860789, 11860789, 11860811, 11860811,
+ + 11860837, 11860837, 11860859, 11860859, 11860867, 11860867,
+ + 11860889, 11860889, 11860897, 11860897, 11860963, 11860963,
+ + 11860969, 11860969, 11860973, 11860973, 11860993, 11860993,
+ + 11861011, 11861011, 11861033, 11861033, 11861071, 11861071 /
+ data (iparam(1,i),iparam(2,i),i= 870, 899) /
+ + 11861081, 11861081, 11861089, 11861089, 11861093, 11861093,
+ + 11861099, 11861099, 11861107, 11861107, 11861131, 11861131,
+ + 11861141, 11861141, 11861159, 11861159, 11861167, 11861167,
+ + 11861191, 11861191, 11861197, 11861197, 11861207, 11861207,
+ + 11861219, 11861219, 11861221, 11861221, 11861231, 11861231,
+ + 11861237, 11861237, 11861273, 11861273, 11861293, 11861293,
+ + 11861299, 11861299, 11861303, 11861303, 11861327, 11861327,
+ + 11861351, 11861351, 11861357, 11861357, 11861363, 11861363,
+ + 11861371, 11861371, 11861401, 11861401, 11861407, 11861407,
+ + 11861411, 11861411, 11861413, 11861413, 11861429, 11861429 /
+ data (iparam(1,i),iparam(2,i),i= 900, 929) /
+ + 11861441, 11861441, 11861467, 11861467, 11861527, 11861527,
+ + 11861539, 11861539, 11861543, 11861543, 11861557, 11861557,
+ + 11861569, 11861569, 11861573, 11861573, 11861579, 11861579,
+ + 11861581, 11861581, 11861599, 11861599, 11861611, 11861611,
+ + 11861617, 11861617, 11861627, 11861627, 11861639, 11861639,
+ + 11861651, 11861651, 11861659, 11861659, 11861671, 11861671,
+ + 11861683, 11861683, 11861687, 11861687, 11861693, 11861693,
+ + 11861701, 11861701, 11861711, 11861711, 11861713, 11861713,
+ + 11861749, 11861749, 11861791, 11861791, 11861803, 11861803,
+ + 11861819, 11861819, 11861827, 11861827, 11861849, 11861849 /
+ data (iparam(1,i),iparam(2,i),i= 930, 959) /
+ + 11861873, 11861873, 11861879, 11861879, 11861887, 11861887,
+ + 11861911, 11861911, 11861917, 11861917, 11861921, 11861921,
+ + 11861923, 11861923, 11861953, 11861953, 11861959, 11861959,
+ + 11861987, 11861987, 11862007, 11862007, 11862013, 11862013,
+ + 11862029, 11862029, 11862031, 11862031, 11862049, 11862049,
+ + 11862077, 11862077, 11862083, 11862083, 11862157, 11862157,
+ + 11862167, 11862167, 11862199, 11862199, 11862203, 11862203,
+ + 11862217, 11862217, 11862223, 11862223, 11862229, 11862229,
+ + 11862233, 11862233, 11862239, 11862239, 11862241, 11862241,
+ + 11862259, 11862259, 11862269, 11862269, 11862271, 11862271 /
+ data (iparam(1,i),iparam(2,i),i= 960, 989) /
+ + 11862293, 11862293, 11862307, 11862307, 11862313, 11862313,
+ + 11862317, 11862317, 11862343, 11862343, 11862353, 11862353,
+ + 11862373, 11862373, 11862391, 11862391, 11862439, 11862439,
+ + 11862469, 11862469, 11862493, 11862493, 11862527, 11862527,
+ + 11862547, 11862547, 11862563, 11862563, 11862569, 11862569,
+ + 11862577, 11862577, 11862581, 11862581, 11862611, 11862611,
+ + 11862623, 11862623, 11862661, 11862661, 11862673, 11862673,
+ + 11862679, 11862679, 11862701, 11862701, 11862703, 11862703,
+ + 11862713, 11862713, 11862761, 11862761, 11862791, 11862791,
+ + 11862803, 11862803, 11862839, 11862839, 11862841, 11862841 /
+ data (iparam(1,i),iparam(2,i),i= 990,1019) /
+ + 11862857, 11862857, 11862869, 11862869, 11862881, 11862881,
+ + 11862911, 11862911, 11862919, 11862919, 11862959, 11862959,
+ + 11862979, 11862979, 11862989, 11862989, 11862997, 11862997,
+ + 11863021, 11863021, 11863031, 11863031, 11863037, 11863037,
+ + 11863039, 11863039, 11863057, 11863057, 11863067, 11863067,
+ + 11863073, 11863073, 11863099, 11863099, 11863109, 11863109,
+ + 11863121, 11863121, 11863123, 11863123, 11863133, 11863133,
+ + 11863151, 11863151, 11863153, 11863153, 11863171, 11863171,
+ + 11863183, 11863183, 11863207, 11863207, 11863213, 11863213,
+ + 11863237, 11863237, 11863249, 11863249, 11863253, 11863253 /
+ data (iparam(1,i),iparam(2,i),i=1020,1021) /
+ + 11863259, 11863259, 11863279, 11863279 /
+ end
parmread.F
pinorm.f
printmat.f
- prng_32.F
+ prng.F
randgens.f
ran.f
readpdb.F
MP.F
newconf.F
parmread.F
- prng_32.F
+ prng.F
readpdb.F
readrtns_csa.F
rmsd.F
--- /dev/null
+ real*8 function prng_next(me)
+ implicit none
+ integer me
+c
+c Calling sequence:
+c <new random number> = prng_next ( <ordinal of generator desired> )
+c <vector of random #s> = vprng ( <ordinal>, <vector>, <length> )
+c
+c This code is based on a sequential algorithm provided by Mal Kalos.
+c This version uses a single 64-bit word to store the initial seeds
+c and additive constants.
+c A 64-bit floating point number is returned.
+c
+c The array "iparam" is full-word aligned, being padded by zeros to
+c let each generator be on a subpage boundary.
+c That is, rows 1 and 2 in a given column of the array are for real,
+c rows 3-16 are bogus.
+c
+c July 12, 1993: double the number of sequences. We should have been
+c using two packets per seed, rather than four
+c October 31, 1993: merge the two arrays of seeds and constants,
+c and switch to 64-bit arithmetic.
+c June 1994: port to RS6K. Internal state is kept as 2 64-bit integers
+c The ishft function is defined only on 32-bit integers, so we will
+c shift numbers by dividing by 2**11 and then adding on 2**53-1.
+c
+c November 1994: ishift now works on 64-bit numbers (though it gives a
+c warning). Thus we go back to using it. John Zollweg also added the
+c vprng() routine to return vectors of real*8 random numbers.
+c
+ real*8 recip53
+ parameter ( recip53 = 2.0D0**(-53) )
+ integer*8 two
+ parameter ( two = 2**11)
+ integer*8 m,ishift
+c parameter ( m = 34522712143931 ) ! 11**13
+c parameter ( ishift = 9007199254740991 ) ! 2**53-1
+
+ integer nmax
+ integer*8 iparam
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ integer*8 next
+
+crc g77 doesn't support integer*8 constants
+ m = dint(34522712143931.0d0)
+ ishift = dint(9007199254740991.0d0)
+
+c RS6K porting note: ishift now takes 64-bit integers , with a warning
+ if ( 0.le.me .and. me.le.nmax ) then
+ next = iparam(1,me)*m + iparam(2,me)
+ iparam(1,me) = next
+ prng_next = recip53 * ishft( next, -11 )
+ else
+ prng_next=-1.0D0
+ endif
+
+ end
+c
+c vprng(me, rn, num) Get a vector of random numbers
+c
+ subroutine vprng(me,rn,num)
+ real*8 recip53, rn(1)
+ parameter ( recip53 = 2.0D0**(-53) )
+ integer*8 m,iparam
+c parameter ( m = 34522712143931 ) ! 11**13
+ integer nmax, num, me
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ integer*8 next
+
+crc g77 doesn't support integer*8 constants
+ m = dint(34522712143931.0d0)
+
+ if ( 0.le.me .and. me.le.nmax ) then
+ do 1 i=1,num
+ next = iparam(1,me)*m + iparam(2,me)
+ iparam(1,me) = next
+ rn(i) = recip53 * ishft( next, -11 )
+ 1 continue
+ else
+ rn(1)=-1.0D0
+ endif
+ return
+ end
+
+c
+c prng_chkpnt Get the current state of a generator
+c
+c Calling sequence:
+c logical prng_chkpnt, status
+c status = prng_chkpnt (me, iseed) where
+c
+c me is the particular generator whose state is being gotten
+c seed is an 4-element integer array where the "l"-values will be saved
+c
+ logical function prng_chkpnt (me, iseed)
+ implicit none
+ integer me
+ integer*8 iseed
+
+ integer nmax
+ integer*8 iparam
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ if (me .lt. 0 .or. me .gt. nmax) then
+ prng_chkpnt=.false.
+ else
+ prng_chkpnt=.true.
+ iseed=iparam(1,me)
+ endif
+ end
+c
+c prng_restart Restart generator from a saved state
+c
+c Calling sequence:
+c logical prng_restart, status
+c status = prng_restart (me, iseed) where
+c
+c me is the particular generator being restarted
+c iseed is a 8-byte integer containing the "l"-values
+c
+ logical function prng_restart (me, iseed)
+ implicit none
+ integer me
+ integer*8 iseed
+
+ integer nmax
+ integer*8 iparam
+ parameter(nmax=1021)
+ common/ksrprng/iparam(2,0:nmax)
+
+ if (me .lt. 0 .or. me .gt. nmax) then
+ prng_restart=.false.
+ return
+ else
+ prng_restart=.true.
+ iparam(1,me)=iseed
+ endif
+ end
+
+ block data prngblk
+ parameter(nmax=1021)
+ integer*8 iparam
+ common/ksrprng/iparam(2,0:nmax)
+ data (iparam(1,i),iparam(2,i),i= 0, 29) /
+ + 11848219, 11848219, 11848237, 11848237, 11848241, 11848241,
+ + 11848247, 11848247, 11848253, 11848253, 11848271, 11848271,
+ + 11848297, 11848297, 11848313, 11848313, 11848339, 11848339,
+ + 11848351, 11848351, 11848357, 11848357, 11848363, 11848363,
+ + 11848367, 11848367, 11848373, 11848373, 11848379, 11848379,
+ + 11848393, 11848393, 11848433, 11848433, 11848451, 11848451,
+ + 11848469, 11848469, 11848477, 11848477, 11848489, 11848489,
+ + 11848493, 11848493, 11848513, 11848513, 11848523, 11848523,
+ + 11848531, 11848531, 11848537, 11848537, 11848553, 11848553,
+ + 11848589, 11848589, 11848591, 11848591, 11848601, 11848601 /
+ data (iparam(1,i),iparam(2,i),i= 30, 59) /
+ + 11848619, 11848619, 11848637, 11848637, 11848663, 11848663,
+ + 11848673, 11848673, 11848679, 11848679, 11848691, 11848691,
+ + 11848699, 11848699, 11848709, 11848709, 11848717, 11848717,
+ + 11848721, 11848721, 11848729, 11848729, 11848741, 11848741,
+ + 11848751, 11848751, 11848757, 11848757, 11848787, 11848787,
+ + 11848801, 11848801, 11848829, 11848829, 11848853, 11848853,
+ + 11848861, 11848861, 11848867, 11848867, 11848873, 11848873,
+ + 11848891, 11848891, 11848909, 11848909, 11848919, 11848919,
+ + 11848931, 11848931, 11848937, 11848937, 11848961, 11848961,
+ + 11848981, 11848981, 11849021, 11849021, 11849039, 11849039 /
+ data (iparam(1,i),iparam(2,i),i= 60, 89) /
+ + 11849053, 11849053, 11849059, 11849059, 11849069, 11849069,
+ + 11849077, 11849077, 11849087, 11849087, 11849093, 11849093,
+ + 11849107, 11849107, 11849111, 11849111, 11849129, 11849129,
+ + 11849137, 11849137, 11849177, 11849177, 11849183, 11849183,
+ + 11849203, 11849203, 11849231, 11849231, 11849237, 11849237,
+ + 11849239, 11849239, 11849249, 11849249, 11849251, 11849251,
+ + 11849269, 11849269, 11849273, 11849273, 11849291, 11849291,
+ + 11849297, 11849297, 11849309, 11849309, 11849339, 11849339,
+ + 11849359, 11849359, 11849363, 11849363, 11849399, 11849399,
+ + 11849401, 11849401, 11849413, 11849413, 11849417, 11849417 /
+ data (iparam(1,i),iparam(2,i),i= 90, 119) /
+ + 11849437, 11849437, 11849443, 11849443, 11849473, 11849473,
+ + 11849491, 11849491, 11849503, 11849503, 11849507, 11849507,
+ + 11849557, 11849557, 11849567, 11849567, 11849569, 11849569,
+ + 11849573, 11849573, 11849587, 11849587, 11849599, 11849599,
+ + 11849633, 11849633, 11849641, 11849641, 11849653, 11849653,
+ + 11849659, 11849659, 11849671, 11849671, 11849683, 11849683,
+ + 11849689, 11849689, 11849693, 11849693, 11849699, 11849699,
+ + 11849701, 11849701, 11849707, 11849707, 11849713, 11849713,
+ + 11849723, 11849723, 11849741, 11849741, 11849743, 11849743,
+ + 11849759, 11849759, 11849767, 11849767, 11849771, 11849771 /
+ data (iparam(1,i),iparam(2,i),i= 120, 149) /
+ + 11849791, 11849791, 11849801, 11849801, 11849809, 11849809,
+ + 11849813, 11849813, 11849869, 11849869, 11849881, 11849881,
+ + 11849891, 11849891, 11849909, 11849909, 11849923, 11849923,
+ + 11849933, 11849933, 11849947, 11849947, 11849987, 11849987,
+ + 11850001, 11850001, 11850011, 11850011, 11850019, 11850019,
+ + 11850023, 11850023, 11850031, 11850031, 11850049, 11850049,
+ + 11850061, 11850061, 11850073, 11850073, 11850077, 11850077,
+ + 11850103, 11850103, 11850109, 11850109, 11850121, 11850121,
+ + 11850127, 11850127, 11850133, 11850133, 11850149, 11850149,
+ + 11850161, 11850161, 11850169, 11850169, 11850191, 11850191 /
+ data (iparam(1,i),iparam(2,i),i= 150, 179) /
+ + 11850233, 11850233, 11850247, 11850247, 11850259, 11850259,
+ + 11850269, 11850269, 11850283, 11850283, 11850301, 11850301,
+ + 11850341, 11850341, 11850347, 11850347, 11850367, 11850367,
+ + 11850373, 11850373, 11850379, 11850379, 11850389, 11850389,
+ + 11850407, 11850407, 11850427, 11850427, 11850437, 11850437,
+ + 11850469, 11850469, 11850481, 11850481, 11850511, 11850511,
+ + 11850529, 11850529, 11850541, 11850541, 11850557, 11850557,
+ + 11850607, 11850607, 11850611, 11850611, 11850667, 11850667,
+ + 11850677, 11850677, 11850679, 11850679, 11850701, 11850701,
+ + 11850731, 11850731, 11850739, 11850739, 11850749, 11850749 /
+ data (iparam(1,i),iparam(2,i),i= 180, 209) /
+ + 11850791, 11850791, 11850803, 11850803, 11850829, 11850829,
+ + 11850833, 11850833, 11850859, 11850859, 11850877, 11850877,
+ + 11850899, 11850899, 11850907, 11850907, 11850913, 11850913,
+ + 11850919, 11850919, 11850931, 11850931, 11850941, 11850941,
+ + 11850947, 11850947, 11850953, 11850953, 11850961, 11850961,
+ + 11850983, 11850983, 11850991, 11850991, 11850997, 11850997,
+ + 11851031, 11851031, 11851033, 11851033, 11851051, 11851051,
+ + 11851061, 11851061, 11851067, 11851067, 11851093, 11851093,
+ + 11851109, 11851109, 11851123, 11851123, 11851127, 11851127,
+ + 11851139, 11851139, 11851157, 11851157, 11851163, 11851163 /
+ data (iparam(1,i),iparam(2,i),i= 210, 239) /
+ + 11851181, 11851181, 11851201, 11851201, 11851219, 11851219,
+ + 11851291, 11851291, 11851303, 11851303, 11851309, 11851309,
+ + 11851313, 11851313, 11851319, 11851319, 11851349, 11851349,
+ + 11851351, 11851351, 11851361, 11851361, 11851373, 11851373,
+ + 11851403, 11851403, 11851409, 11851409, 11851423, 11851423,
+ + 11851447, 11851447, 11851451, 11851451, 11851481, 11851481,
+ + 11851493, 11851493, 11851519, 11851519, 11851523, 11851523,
+ + 11851529, 11851529, 11851547, 11851547, 11851549, 11851549,
+ + 11851559, 11851559, 11851577, 11851577, 11851589, 11851589,
+ + 11851591, 11851591, 11851597, 11851597, 11851603, 11851603 /
+ data (iparam(1,i),iparam(2,i),i= 240, 269) /
+ + 11851607, 11851607, 11851613, 11851613, 11851621, 11851621,
+ + 11851627, 11851627, 11851639, 11851639, 11851673, 11851673,
+ + 11851681, 11851681, 11851727, 11851727, 11851753, 11851753,
+ + 11851759, 11851759, 11851787, 11851787, 11851793, 11851793,
+ + 11851799, 11851799, 11851813, 11851813, 11851841, 11851841,
+ + 11851859, 11851859, 11851867, 11851867, 11851891, 11851891,
+ + 11851909, 11851909, 11851919, 11851919, 11851927, 11851927,
+ + 11851933, 11851933, 11851949, 11851949, 11851967, 11851967,
+ + 11851997, 11851997, 11852017, 11852017, 11852051, 11852051,
+ + 11852053, 11852053, 11852059, 11852059, 11852083, 11852083 /
+ data (iparam(1,i),iparam(2,i),i= 270, 299) /
+ + 11852089, 11852089, 11852129, 11852129, 11852147, 11852147,
+ + 11852149, 11852149, 11852161, 11852161, 11852171, 11852171,
+ + 11852177, 11852177, 11852209, 11852209, 11852221, 11852221,
+ + 11852237, 11852237, 11852251, 11852251, 11852263, 11852263,
+ + 11852273, 11852273, 11852279, 11852279, 11852287, 11852287,
+ + 11852293, 11852293, 11852297, 11852297, 11852303, 11852303,
+ + 11852311, 11852311, 11852327, 11852327, 11852339, 11852339,
+ + 11852341, 11852341, 11852359, 11852359, 11852369, 11852369,
+ + 11852437, 11852437, 11852453, 11852453, 11852459, 11852459,
+ + 11852473, 11852473, 11852513, 11852513, 11852531, 11852531 /
+ data (iparam(1,i),iparam(2,i),i= 300, 329) /
+ + 11852537, 11852537, 11852539, 11852539, 11852557, 11852557,
+ + 11852573, 11852573, 11852579, 11852579, 11852591, 11852591,
+ + 11852609, 11852609, 11852611, 11852611, 11852623, 11852623,
+ + 11852641, 11852641, 11852647, 11852647, 11852657, 11852657,
+ + 11852663, 11852663, 11852717, 11852717, 11852719, 11852719,
+ + 11852741, 11852741, 11852759, 11852759, 11852767, 11852767,
+ + 11852773, 11852773, 11852803, 11852803, 11852807, 11852807,
+ + 11852809, 11852809, 11852831, 11852831, 11852833, 11852833,
+ + 11852837, 11852837, 11852857, 11852857, 11852873, 11852873,
+ + 11852879, 11852879, 11852891, 11852891, 11852917, 11852917 /
+ data (iparam(1,i),iparam(2,i),i= 330, 359) /
+ + 11852921, 11852921, 11852957, 11852957, 11852959, 11852959,
+ + 11852969, 11852969, 11852983, 11852983, 11852989, 11852989,
+ + 11853001, 11853001, 11853013, 11853013, 11853019, 11853019,
+ + 11853031, 11853031, 11853089, 11853089, 11853133, 11853133,
+ + 11853157, 11853157, 11853161, 11853161, 11853181, 11853181,
+ + 11853203, 11853203, 11853217, 11853217, 11853221, 11853221,
+ + 11853227, 11853227, 11853241, 11853241, 11853307, 11853307,
+ + 11853319, 11853319, 11853323, 11853323, 11853329, 11853329,
+ + 11853367, 11853367, 11853383, 11853383, 11853419, 11853419,
+ + 11853421, 11853421, 11853427, 11853427, 11853449, 11853449 /
+ data (iparam(1,i),iparam(2,i),i= 360, 389) /
+ + 11853451, 11853451, 11853463, 11853463, 11853529, 11853529,
+ + 11853557, 11853557, 11853571, 11853571, 11853601, 11853601,
+ + 11853613, 11853613, 11853617, 11853617, 11853629, 11853629,
+ + 11853649, 11853649, 11853659, 11853659, 11853679, 11853679,
+ + 11853689, 11853689, 11853719, 11853719, 11853731, 11853731,
+ + 11853757, 11853757, 11853761, 11853761, 11853773, 11853773,
+ + 11853791, 11853791, 11853817, 11853817, 11853839, 11853839,
+ + 11853847, 11853847, 11853857, 11853857, 11853869, 11853869,
+ + 11853883, 11853883, 11853887, 11853887, 11853889, 11853889,
+ + 11853893, 11853893, 11853899, 11853899, 11853911, 11853911 /
+ data (iparam(1,i),iparam(2,i),i= 390, 419) /
+ + 11853931, 11853931, 11853943, 11853943, 11853979, 11853979,
+ + 11853991, 11853991, 11854001, 11854001, 11854009, 11854009,
+ + 11854019, 11854019, 11854057, 11854057, 11854061, 11854061,
+ + 11854147, 11854147, 11854159, 11854159, 11854163, 11854163,
+ + 11854169, 11854169, 11854211, 11854211, 11854247, 11854247,
+ + 11854261, 11854261, 11854267, 11854267, 11854279, 11854279,
+ + 11854303, 11854303, 11854327, 11854327, 11854331, 11854331,
+ + 11854333, 11854333, 11854363, 11854363, 11854379, 11854379,
+ + 11854399, 11854399, 11854411, 11854411, 11854429, 11854429,
+ + 11854433, 11854433, 11854439, 11854439, 11854441, 11854441 /
+ data (iparam(1,i),iparam(2,i),i= 420, 449) /
+ + 11854463, 11854463, 11854477, 11854477, 11854489, 11854489,
+ + 11854517, 11854517, 11854519, 11854519, 11854523, 11854523,
+ + 11854529, 11854529, 11854567, 11854567, 11854571, 11854571,
+ + 11854573, 11854573, 11854603, 11854603, 11854607, 11854607,
+ + 11854681, 11854681, 11854691, 11854691, 11854709, 11854709,
+ + 11854723, 11854723, 11854757, 11854757, 11854783, 11854783,
+ + 11854793, 11854793, 11854813, 11854813, 11854847, 11854847,
+ + 11854853, 11854853, 11854873, 11854873, 11854877, 11854877,
+ + 11854883, 11854883, 11854891, 11854891, 11854897, 11854897,
+ + 11854901, 11854901, 11854919, 11854919, 11854937, 11854937 /
+ data (iparam(1,i),iparam(2,i),i= 450, 479) /
+ + 11854961, 11854961, 11854963, 11854963, 11854979, 11854979,
+ + 11855003, 11855003, 11855017, 11855017, 11855023, 11855023,
+ + 11855029, 11855029, 11855033, 11855033, 11855111, 11855111,
+ + 11855141, 11855141, 11855147, 11855147, 11855149, 11855149,
+ + 11855159, 11855159, 11855177, 11855177, 11855203, 11855203,
+ + 11855213, 11855213, 11855219, 11855219, 11855231, 11855231,
+ + 11855267, 11855267, 11855269, 11855269, 11855303, 11855303,
+ + 11855309, 11855309, 11855321, 11855321, 11855329, 11855329,
+ + 11855339, 11855339, 11855351, 11855351, 11855353, 11855353,
+ + 11855357, 11855357, 11855359, 11855359, 11855381, 11855381 /
+ data (iparam(1,i),iparam(2,i),i= 480, 509) /
+ + 11855383, 11855383, 11855387, 11855387, 11855399, 11855399,
+ + 11855407, 11855407, 11855413, 11855413, 11855489, 11855489,
+ + 11855491, 11855491, 11855507, 11855507, 11855521, 11855521,
+ + 11855531, 11855531, 11855549, 11855549, 11855551, 11855551,
+ + 11855567, 11855567, 11855581, 11855581, 11855587, 11855587,
+ + 11855593, 11855593, 11855633, 11855633, 11855653, 11855653,
+ + 11855663, 11855663, 11855687, 11855687, 11855689, 11855689,
+ + 11855699, 11855699, 11855713, 11855713, 11855731, 11855731,
+ + 11855737, 11855737, 11855743, 11855743, 11855747, 11855747,
+ + 11855759, 11855759, 11855773, 11855773, 11855801, 11855801 /
+ data (iparam(1,i),iparam(2,i),i= 510, 539) /
+ + 11855807, 11855807, 11855813, 11855813, 11855827, 11855827,
+ + 11855839, 11855839, 11855869, 11855869, 11855881, 11855881,
+ + 11855903, 11855903, 11855911, 11855911, 11855933, 11855933,
+ + 11855959, 11855959, 11855989, 11855989, 11855993, 11855993,
+ + 11855999, 11855999, 11856001, 11856001, 11856023, 11856023,
+ + 11856049, 11856049, 11856071, 11856071, 11856101, 11856101,
+ + 11856107, 11856107, 11856113, 11856113, 11856139, 11856139,
+ + 11856151, 11856151, 11856161, 11856161, 11856179, 11856179,
+ + 11856193, 11856193, 11856199, 11856199, 11856223, 11856223,
+ + 11856239, 11856239, 11856263, 11856263, 11856269, 11856269 /
+ data (iparam(1,i),iparam(2,i),i= 540, 569) /
+ + 11856281, 11856281, 11856287, 11856287, 11856307, 11856307,
+ + 11856311, 11856311, 11856329, 11856329, 11856343, 11856343,
+ + 11856359, 11856359, 11856371, 11856371, 11856373, 11856373,
+ + 11856409, 11856409, 11856419, 11856419, 11856461, 11856461,
+ + 11856469, 11856469, 11856473, 11856473, 11856479, 11856479,
+ + 11856511, 11856511, 11856517, 11856517, 11856541, 11856541,
+ + 11856547, 11856547, 11856553, 11856553, 11856583, 11856583,
+ + 11856629, 11856629, 11856641, 11856641, 11856653, 11856653,
+ + 11856659, 11856659, 11856673, 11856673, 11856697, 11856697,
+ + 11856709, 11856709, 11856727, 11856727, 11856731, 11856731 /
+ data (iparam(1,i),iparam(2,i),i= 570, 599) /
+ + 11856763, 11856763, 11856809, 11856809, 11856811, 11856811,
+ + 11856821, 11856821, 11856841, 11856841, 11856857, 11856857,
+ + 11856877, 11856877, 11856883, 11856883, 11856899, 11856899,
+ + 11856919, 11856919, 11856947, 11856947, 11856953, 11856953,
+ + 11856979, 11856979, 11857003, 11857003, 11857033, 11857033,
+ + 11857037, 11857037, 11857039, 11857039, 11857049, 11857049,
+ + 11857061, 11857061, 11857067, 11857067, 11857073, 11857073,
+ + 11857081, 11857081, 11857091, 11857091, 11857093, 11857093,
+ + 11857099, 11857099, 11857123, 11857123, 11857127, 11857127,
+ + 11857147, 11857147, 11857151, 11857151, 11857193, 11857193 /
+ data (iparam(1,i),iparam(2,i),i= 600, 629) /
+ + 11857217, 11857217, 11857229, 11857229, 11857243, 11857243,
+ + 11857249, 11857249, 11857267, 11857267, 11857277, 11857277,
+ + 11857291, 11857291, 11857303, 11857303, 11857309, 11857309,
+ + 11857327, 11857327, 11857331, 11857331, 11857333, 11857333,
+ + 11857361, 11857361, 11857367, 11857367, 11857369, 11857369,
+ + 11857393, 11857393, 11857399, 11857399, 11857409, 11857409,
+ + 11857421, 11857421, 11857423, 11857423, 11857451, 11857451,
+ + 11857453, 11857453, 11857457, 11857457, 11857477, 11857477,
+ + 11857481, 11857481, 11857493, 11857493, 11857499, 11857499,
+ + 11857519, 11857519, 11857523, 11857523, 11857529, 11857529 /
+ data (iparam(1,i),iparam(2,i),i= 630, 659) /
+ + 11857543, 11857543, 11857561, 11857561, 11857589, 11857589,
+ + 11857591, 11857591, 11857613, 11857613, 11857621, 11857621,
+ + 11857661, 11857661, 11857667, 11857667, 11857693, 11857693,
+ + 11857697, 11857697, 11857709, 11857709, 11857711, 11857711,
+ + 11857751, 11857751, 11857753, 11857753, 11857759, 11857759,
+ + 11857763, 11857763, 11857777, 11857777, 11857787, 11857787,
+ + 11857793, 11857793, 11857801, 11857801, 11857817, 11857817,
+ + 11857819, 11857819, 11857831, 11857831, 11857837, 11857837,
+ + 11857873, 11857873, 11857877, 11857877, 11857883, 11857883,
+ + 11857889, 11857889, 11857907, 11857907, 11857913, 11857913 /
+ data (iparam(1,i),iparam(2,i),i= 660, 689) /
+ + 11857931, 11857931, 11857969, 11857969, 11857991, 11857991,
+ + 11857999, 11857999, 11858009, 11858009, 11858017, 11858017,
+ + 11858023, 11858023, 11858029, 11858029, 11858039, 11858039,
+ + 11858051, 11858051, 11858057, 11858057, 11858059, 11858059,
+ + 11858101, 11858101, 11858111, 11858111, 11858131, 11858131,
+ + 11858149, 11858149, 11858159, 11858159, 11858177, 11858177,
+ + 11858191, 11858191, 11858201, 11858201, 11858227, 11858227,
+ + 11858243, 11858243, 11858267, 11858267, 11858269, 11858269,
+ + 11858279, 11858279, 11858281, 11858281, 11858291, 11858291,
+ + 11858311, 11858311, 11858323, 11858323, 11858359, 11858359 /
+ data (iparam(1,i),iparam(2,i),i= 690, 719) /
+ + 11858377, 11858377, 11858381, 11858381, 11858387, 11858387,
+ + 11858423, 11858423, 11858443, 11858443, 11858447, 11858447,
+ + 11858479, 11858479, 11858533, 11858533, 11858543, 11858543,
+ + 11858551, 11858551, 11858557, 11858557, 11858569, 11858569,
+ + 11858573, 11858573, 11858579, 11858579, 11858597, 11858597,
+ + 11858599, 11858599, 11858629, 11858629, 11858657, 11858657,
+ + 11858659, 11858659, 11858683, 11858683, 11858701, 11858701,
+ + 11858719, 11858719, 11858723, 11858723, 11858729, 11858729,
+ + 11858747, 11858747, 11858779, 11858779, 11858783, 11858783,
+ + 11858801, 11858801, 11858807, 11858807, 11858813, 11858813 /
+ data (iparam(1,i),iparam(2,i),i= 720, 749) /
+ + 11858839, 11858839, 11858851, 11858851, 11858893, 11858893,
+ + 11858897, 11858897, 11858921, 11858921, 11858947, 11858947,
+ + 11858953, 11858953, 11858969, 11858969, 11858971, 11858971,
+ + 11858989, 11858989, 11859017, 11859017, 11859031, 11859031,
+ + 11859049, 11859049, 11859061, 11859061, 11859073, 11859073,
+ + 11859077, 11859077, 11859079, 11859079, 11859083, 11859083,
+ + 11859101, 11859101, 11859109, 11859109, 11859137, 11859137,
+ + 11859139, 11859139, 11859151, 11859151, 11859157, 11859157,
+ + 11859163, 11859163, 11859167, 11859167, 11859179, 11859179,
+ + 11859187, 11859187, 11859229, 11859229, 11859233, 11859233 /
+ data (iparam(1,i),iparam(2,i),i= 750, 779) /
+ + 11859241, 11859241, 11859247, 11859247, 11859269, 11859269,
+ + 11859293, 11859293, 11859307, 11859307, 11859311, 11859311,
+ + 11859349, 11859349, 11859359, 11859359, 11859371, 11859371,
+ + 11859377, 11859377, 11859383, 11859383, 11859427, 11859427,
+ + 11859433, 11859433, 11859451, 11859451, 11859457, 11859457,
+ + 11859461, 11859461, 11859473, 11859473, 11859481, 11859481,
+ + 11859487, 11859487, 11859493, 11859493, 11859503, 11859503,
+ + 11859509, 11859509, 11859539, 11859539, 11859541, 11859541,
+ + 11859563, 11859563, 11859569, 11859569, 11859571, 11859571,
+ + 11859583, 11859583, 11859599, 11859599, 11859611, 11859611 /
+ data (iparam(1,i),iparam(2,i),i= 780, 809) /
+ + 11859643, 11859643, 11859707, 11859707, 11859713, 11859713,
+ + 11859719, 11859719, 11859739, 11859739, 11859751, 11859751,
+ + 11859791, 11859791, 11859817, 11859817, 11859821, 11859821,
+ + 11859833, 11859833, 11859847, 11859847, 11859853, 11859853,
+ + 11859877, 11859877, 11859889, 11859889, 11859893, 11859893,
+ + 11859901, 11859901, 11859907, 11859907, 11859917, 11859917,
+ + 11859923, 11859923, 11859929, 11859929, 11859961, 11859961,
+ + 11859979, 11859979, 11859989, 11859989, 11859997, 11859997,
+ + 11860021, 11860021, 11860031, 11860031, 11860039, 11860039,
+ + 11860049, 11860049, 11860081, 11860081, 11860087, 11860087 /
+ data (iparam(1,i),iparam(2,i),i= 810, 839) /
+ + 11860097, 11860097, 11860103, 11860103, 11860109, 11860109,
+ + 11860117, 11860117, 11860133, 11860133, 11860151, 11860151,
+ + 11860171, 11860171, 11860207, 11860207, 11860223, 11860223,
+ + 11860231, 11860231, 11860243, 11860243, 11860267, 11860267,
+ + 11860301, 11860301, 11860307, 11860307, 11860327, 11860327,
+ + 11860379, 11860379, 11860397, 11860397, 11860411, 11860411,
+ + 11860469, 11860469, 11860477, 11860477, 11860483, 11860483,
+ + 11860487, 11860487, 11860489, 11860489, 11860493, 11860493,
+ + 11860517, 11860517, 11860547, 11860547, 11860567, 11860567,
+ + 11860573, 11860573, 11860613, 11860613, 11860619, 11860619 /
+ data (iparam(1,i),iparam(2,i),i= 840, 869) /
+ + 11860627, 11860627, 11860637, 11860637, 11860643, 11860643,
+ + 11860649, 11860649, 11860661, 11860661, 11860669, 11860669,
+ + 11860687, 11860687, 11860691, 11860691, 11860697, 11860697,
+ + 11860699, 11860699, 11860703, 11860703, 11860727, 11860727,
+ + 11860741, 11860741, 11860753, 11860753, 11860777, 11860777,
+ + 11860787, 11860787, 11860789, 11860789, 11860811, 11860811,
+ + 11860837, 11860837, 11860859, 11860859, 11860867, 11860867,
+ + 11860889, 11860889, 11860897, 11860897, 11860963, 11860963,
+ + 11860969, 11860969, 11860973, 11860973, 11860993, 11860993,
+ + 11861011, 11861011, 11861033, 11861033, 11861071, 11861071 /
+ data (iparam(1,i),iparam(2,i),i= 870, 899) /
+ + 11861081, 11861081, 11861089, 11861089, 11861093, 11861093,
+ + 11861099, 11861099, 11861107, 11861107, 11861131, 11861131,
+ + 11861141, 11861141, 11861159, 11861159, 11861167, 11861167,
+ + 11861191, 11861191, 11861197, 11861197, 11861207, 11861207,
+ + 11861219, 11861219, 11861221, 11861221, 11861231, 11861231,
+ + 11861237, 11861237, 11861273, 11861273, 11861293, 11861293,
+ + 11861299, 11861299, 11861303, 11861303, 11861327, 11861327,
+ + 11861351, 11861351, 11861357, 11861357, 11861363, 11861363,
+ + 11861371, 11861371, 11861401, 11861401, 11861407, 11861407,
+ + 11861411, 11861411, 11861413, 11861413, 11861429, 11861429 /
+ data (iparam(1,i),iparam(2,i),i= 900, 929) /
+ + 11861441, 11861441, 11861467, 11861467, 11861527, 11861527,
+ + 11861539, 11861539, 11861543, 11861543, 11861557, 11861557,
+ + 11861569, 11861569, 11861573, 11861573, 11861579, 11861579,
+ + 11861581, 11861581, 11861599, 11861599, 11861611, 11861611,
+ + 11861617, 11861617, 11861627, 11861627, 11861639, 11861639,
+ + 11861651, 11861651, 11861659, 11861659, 11861671, 11861671,
+ + 11861683, 11861683, 11861687, 11861687, 11861693, 11861693,
+ + 11861701, 11861701, 11861711, 11861711, 11861713, 11861713,
+ + 11861749, 11861749, 11861791, 11861791, 11861803, 11861803,
+ + 11861819, 11861819, 11861827, 11861827, 11861849, 11861849 /
+ data (iparam(1,i),iparam(2,i),i= 930, 959) /
+ + 11861873, 11861873, 11861879, 11861879, 11861887, 11861887,
+ + 11861911, 11861911, 11861917, 11861917, 11861921, 11861921,
+ + 11861923, 11861923, 11861953, 11861953, 11861959, 11861959,
+ + 11861987, 11861987, 11862007, 11862007, 11862013, 11862013,
+ + 11862029, 11862029, 11862031, 11862031, 11862049, 11862049,
+ + 11862077, 11862077, 11862083, 11862083, 11862157, 11862157,
+ + 11862167, 11862167, 11862199, 11862199, 11862203, 11862203,
+ + 11862217, 11862217, 11862223, 11862223, 11862229, 11862229,
+ + 11862233, 11862233, 11862239, 11862239, 11862241, 11862241,
+ + 11862259, 11862259, 11862269, 11862269, 11862271, 11862271 /
+ data (iparam(1,i),iparam(2,i),i= 960, 989) /
+ + 11862293, 11862293, 11862307, 11862307, 11862313, 11862313,
+ + 11862317, 11862317, 11862343, 11862343, 11862353, 11862353,
+ + 11862373, 11862373, 11862391, 11862391, 11862439, 11862439,
+ + 11862469, 11862469, 11862493, 11862493, 11862527, 11862527,
+ + 11862547, 11862547, 11862563, 11862563, 11862569, 11862569,
+ + 11862577, 11862577, 11862581, 11862581, 11862611, 11862611,
+ + 11862623, 11862623, 11862661, 11862661, 11862673, 11862673,
+ + 11862679, 11862679, 11862701, 11862701, 11862703, 11862703,
+ + 11862713, 11862713, 11862761, 11862761, 11862791, 11862791,
+ + 11862803, 11862803, 11862839, 11862839, 11862841, 11862841 /
+ data (iparam(1,i),iparam(2,i),i= 990,1019) /
+ + 11862857, 11862857, 11862869, 11862869, 11862881, 11862881,
+ + 11862911, 11862911, 11862919, 11862919, 11862959, 11862959,
+ + 11862979, 11862979, 11862989, 11862989, 11862997, 11862997,
+ + 11863021, 11863021, 11863031, 11863031, 11863037, 11863037,
+ + 11863039, 11863039, 11863057, 11863057, 11863067, 11863067,
+ + 11863073, 11863073, 11863099, 11863099, 11863109, 11863109,
+ + 11863121, 11863121, 11863123, 11863123, 11863133, 11863133,
+ + 11863151, 11863151, 11863153, 11863153, 11863171, 11863171,
+ + 11863183, 11863183, 11863207, 11863207, 11863213, 11863213,
+ + 11863237, 11863237, 11863249, 11863249, 11863253, 11863253 /
+ data (iparam(1,i),iparam(2,i),i=1020,1021) /
+ + 11863259, 11863259, 11863279, 11863279 /
+ end
& g_corr5_loc(maxvar),g_corr6_loc(maxvar),gsccorc(3,maxres),
& gsccorx(3,maxres),gsccor_loc(maxres),dtheta(3,2,maxres),
& gscloc(3,maxres),gsclocx(3,maxres),
- & dphi(3,3,maxres),dalpha(3,3,maxres),domega(3,3,maxres),nfl,icg
+ & dphi(3,3,maxres),dalpha(3,3,maxres),domega(3,3,maxres),nfl,icg,
+ & gdfad(3,maxres),gdfat(3,maxres),gdfan(3,maxres),gdfab(3,maxres)
+
double precision derx,derx_turn
common /deriv_loc/ derx(3,5,2),derx_turn(3,5,2)
double precision dXX_C1tab(3,maxres),dYY_C1tab(3,maxres),
--- /dev/null
+C =======
+C COMMON.DFA
+C =======
+C 2010/12/20 By Juyong Lee
+C
+c parameter
+C [ 8 * ( Nres - 8 ) ] distance restraints
+C [ 2 * ( Nres - 8 ) ] angle restraints
+C [ Nres ] neighbor restraints
+C Total : ~ 11 * Nres restraints
+C
+C
+ INTEGER IDFAMAX,IDFAMX2,IDFACMD,IDMAXMIN, MAXN
+ PARAMETER(IDFAMAX=4000,IDFAMX2=1000,IDFACMD=500,IDMAXMIN=500)
+ PARAMETER(MAXN=4)
+ real*8 wwdist,wwangle,wwnei
+ parameter(wwdist=1.0d0,wwangle=1.0d0,wwnei=1.0d0)
+
+C IDFAMAX - maximum number of DFA restraint including distance, angle and
+C number of neighbors ( Max of assign statement )
+C IDFAMX2 - maximum number of atoms which are targets of restraints
+C IDFACMD - maximum number of 'DFA' command call
+C IDMAXMIN - Maximum number of minima of dist, angle and neighbor info. from fragments
+C MAXN - Maximum Number of shell, currently 4
+C MAXRES - Maximum number of CAs
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
+C INTEGER
+C DFANUM - Number of ALL DFA restrants
+c IDFA[DIS, PHI, THE, NEI] - NUMBER of restraints
+c IDISNUM - number of minima for a distance restraint
+c IPHINUM - number of minima for a phi angle restraint
+c ITHENUM - number of minima for a theta angle restraint
+c INEINUM - number of minima for a number of neighbors restraint
+
+c IDISLIS - atom number of two atoms for distance restraint
+c IPHILIS - atom numbers of four atoms for angle restraint
+c ITHELIS - atom numbers of four atoms for angle restraint
+c INEILIS - atom number of center of neighbor calculation
+c JNEILIS - atom number of target of neighboring calculation
+c JNEINUM - number of target atoms of neighboring term
+C KSHELL - SHELL number
+
+C ishiftca - index shift for CA atoms in UNRES (1 if the 1st aa != GLY)
+C ilastca - index of the last CA atom in UNRES (nres-1 if last aa != GLY)
+
+C old only for CHARMM
+C STOAGDF - Store assign information ( How many assign within one command )
+C NMAP - mapping between dfanum and ndis, nphi, nthe, nnei
+
+ INTEGER IDFADIS,IDFAPHI,IDFATHE,IDFANEI,
+ & IDISLIS,IPHILIS,ITHELIS,INEILIS,
+ & IDISNUM,IPHINUM,ITHENUM,INEINUM,
+ & FNEI,
+ & NCA,ICAIDX,
+ & STOAGDF, NMAP, IDFACAT, KDISNUM, KSHELL
+ & ishiftca,ilastca
+ COMMON /IDFA/ DFACMD, DFANUM,
+ & IDFADIS, IDFAPHI, IDFANEI, IDFATHE,
+ & IDISNUM(IDFAMAX), IPHINUM(IDFAMAX),
+ & ITHENUM(IDFAMAX), INEINUM(IDFAMAX),
+ & FNEI(IDFAMAX,IDMAXMIN), IDISLIS(2,IDFAMAX),
+ & IPHILIS(5,IDFAMAX), ITHELIS(5,IDFAMAX),
+ & INEILIS(IDFAMAX),
+ & KSHELL(IDFAMAX),
+ & IDFACAT(IDFACMD),
+ & KDISNUM(IDFAMAX),
+ & NCA, ICAIDX(MAXRES)
+ COMMON /IDFA2/ ishiftca,ilastca
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+C
+C REAL VARIABLES
+C
+c SCC[DIST, PHI, THE] - weight of each calculations
+c FDIST - distance minima
+C FPHI - phi minima
+c FTHE - theta minima
+C DFAEXP : calculate expential function in advance
+C
+ REAL*8 SCCDIST, SCCPHI, SCCTHE, SCCNEI, FDIST, FPHI1, FPHI2,
+ & FTHE1, FTHE2,
+ & DIS_INC, PHI_INC, THE_INC, NEI_INC, BETA_INC,
+ & WSHET, EDFABET,
+ & CK, SCK
+c & ,DFAEXP
+
+ COMMON /RDFA/ SCCDIST(IDFAMAX,IDMAXMIN),FDIST(IDFAMAX,IDMAXMIN),
+ & SCCPHI(IDFAMAX,IDMAXMIN), SCCTHE(IDFAMAX,IDMAXMIN),
+ & SCCNEI(IDFAMAX,IDMAXMIN),
+ & FPHI1(IDFAMAX,IDMAXMIN), FPHI2(IDFAMAX,IDMAXMIN),
+ & FTHE1(IDFAMAX,IDMAXMIN), FTHE2(IDFAMAX,IDMAXMIN),
+ & DIS_INC, PHI_INC, THE_INC, NEI_INC, BETA_INC,
+ & WSHET(MAXRES,MAXRES), EDFABET,
+ & CK(4),SCK(4),S1(4),S2(4)
+c & ,DFAEXP(15001),
+
+ DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/
+ DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/
+ DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/
+ DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/
common /ffield/ wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,
& wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
& wturn6,wvdwpp,wsct,weights(n_ene),temp0,
+ & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta,
& scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp,
& rescale_mode
common /potentials/ potname(5)
parameter (max_pool=10)
C Number of energy components
integer n_ene,n_ene2
- parameter (n_ene=23,n_ene2=2*n_ene)
+ parameter (n_ene=27,n_ene2=2*n_ene)
C Number of threads in deformation
integer max_thread,max_thread2
parameter (max_thread=4,max_thread2=2*max_thread)
--- /dev/null
+ subroutine init_dfa_vars
+
+ include 'DIMENSIONS'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DFA'
+
+ integer ii
+
+C Number of restraints
+ idisnum = 0
+ iphinum = 0
+ ithenum = 0
+ ineinum = 0
+
+ idislis = 0
+ iphilis = 0
+ ithelis = 0
+ ineilis = 0
+ jneilis = 0
+ jneinum = 0
+ kshell = 0
+ fnei = 0
+C For beta
+ nca = 0
+ icaidx = 0
+
+C real variables
+CC WEIGHTS for each min
+ sccdist = 0.0d0
+ fdist = 0.0d0
+ sccphi = 0.0d0
+ sccthe = 0.0d0
+ sccnei = 0.0d0
+ fphi1 = 0.0d0
+ fphi2 = 0.0d0
+ fthe1 = 0.0d0
+ fthe2 = 0.0d0
+C energies
+ edfatot = 0.0d0
+ edfadis = 0.0d0
+ edfaphi = 0.0d0
+ edfathe = 0.0d0
+ edfanei = 0.0d0
+ edfabet = 0.0d0
+C weights for each E term
+C these should be identical with
+ dis_inc = 0.0d0
+ phi_inc = 0.0d0
+ the_inc = 0.0d0
+ nei_inc = 0.0d0
+ beta_inc = 0.0d0
+ wshet = 0.0d0
+C precalculate exp table!
+c dfaexp = 0.0d0
+c do ii = 1, 15001
+c dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
+c end do
+
+ ishiftca=nnt-1
+ ilastca=nct
+
+ print *,'ishiftca=',ishiftca,'ilastca=',ilastca
+
+ return
+ end
+
+
+ subroutine read_dfa_info
+C
+C read fragment informations
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DFA'
+
+
+C NOTE THAT FILENAMES are FIXED, CURRENTLY!!
+C THIS SHOULD BE MODIFIED!!
+
+ character*320 buffer
+ integer iodfa
+ parameter(iodfa=89)
+
+ integer i, j, nval
+ integer ica1, ica2,ica3,ica4,ica5
+ integer ishell, inca, itmp,iitmp
+ double precision wtmp
+C
+C READ DISTANCE
+C
+ open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
+ goto 34
+ 33 write(iout,'(a)') 'Error opening dist_dfa.dat file'
+ stop
+ 34 continue
+ write(iout,'(a)') 'dist_dfa.dat is opened!'
+C read title
+ read(iodfa, '(a)') buffer
+C read number of restraints
+ read(iodfa, *) IDFADIS
+ read(iodfa, *) dis_inc
+ do i=1, idfadis
+ read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
+
+ idisnum(i)=nval
+ idislis(1,i)=ica1
+ idislis(2,i)=ica2
+
+ do j=1, nval
+ read(iodfa,*) tmp
+ fdist(i,j) = tmp
+ enddo
+
+ do j=1, nval
+ read(iodfa,*) tmp
+ sccdist(i,j) = tmp
+ enddo
+
+ enddo
+ close(iodfa)
+
+C READ ANGLE RESTRAINTS
+C PHI RESTRAINTS
+ open(iodfa, file='phi_dfa.dat',status='old',err=35)
+ goto 36
+ 35 write(iout,'(a)') 'Error opening dist_dfa.dat file'
+ stop
+
+ 36 continue
+ write(iout,'(a)') 'phi_dfa.dat is opened!'
+
+C READ TITLE
+ read(iodfa, '(a)') buffer
+C READ NUMBER OF RESTRAINTS
+ READ(iodfa, *) IDFAPHI
+ read(iodfa,*) phi_inc
+ do i=1, idfaphi
+ read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
+
+ iphinum(i)=nval
+
+ iphilis(1,i)=ica1
+ iphilis(2,i)=ica2
+ iphilis(3,i)=ica3
+ iphilis(4,i)=ica4
+ iphilis(5,i)=ica5
+
+ do j=1, nval
+ read(iodfa,*) tmp1,tmp2
+ fphi1(i,j) = tmp1
+ fphi2(i,j) = tmp2
+ enddo
+
+ do j=1, nval
+ read(iodfa,*) tmp
+ sccphi(i,j) = tmp
+ enddo
+
+ enddo
+ close(iodfa)
+
+C THETA RESTRAINTS
+ open(iodfa, file='theta_dfa.dat',status='old',err=41)
+ goto 42
+ 41 write(iout,'(a)') 'Error opening dist_dfa.dat file'
+ stop
+ 42 continue
+ write(iout,'(a)') 'theta_dfa.dat is opened!'
+C READ TITLE
+ read(iodfa, '(a)') buffer
+C READ NUMBER OF RESTRAINTS
+ READ(iodfa, *) IDFATHE
+ read(iodfa,*) the_inc
+
+ do i=1, idfathe
+ read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
+
+ ithenum(i)=nval
+
+ ithelis(1,i)=ica1
+ ithelis(2,i)=ica2
+ ithelis(3,i)=ica3
+ ithelis(4,i)=ica4
+ ithelis(5,i)=ica5
+
+ do j=1, nval
+ read(iodfa,*) tmp1,tmp2
+ fthe1(i,j) = tmp1
+ fthe2(i,j) = tmp2
+ enddo
+
+ do j=1, nval
+ read(iodfa,*) tmp
+ sccthe(i,j) = tmp
+ enddo
+
+ enddo
+ close(iodfa)
+C END of READING ANGLE RESTRAINT!
+
+C NUMBER OF NEIGHBOR CAs
+ open(iodfa,file='nei_dfa.dat',status='old',err=37)
+ goto 38
+ 37 write(iout,'(a)') 'Error opening nei_dfa.dat file'
+ stop
+ 38 continue
+ write(iout,'(a)') 'nei_dfa.dat is opened!'
+C READ TITLE
+ read(iodfa, '(a)') buffer
+C READ NUMBER OF RESTRAINTS
+ READ(iodfa, *) idfanei
+ read(iodfa,*) nei_inc
+
+ do i=1, idfanei
+ read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval
+
+ ineilis(i)=ica1
+ kshell(i)=ishell
+ ineinum(i)=nval
+
+ do j=1, nval
+ read(iodfa,*) inca
+ fnei(i,j) = inca
+C write(*,*) 'READ NEI:',i,j,fnei(i,j)
+ enddo
+
+ do j=1, nval
+ read(iodfa,*) tmp
+ sccnei(i,j) = tmp
+ enddo
+
+ enddo
+ close(iodfa)
+C END OF NEIGHBORING CA
+
+C READ BETA RESTRAINT
+ open(iodfa, file='beta_dfa.dat',status='old',err=39)
+ goto 40
+ 39 write(iout,'(a)') 'Error opening beta_dfa.dat file'
+ stop
+ 40 continue
+ write(iout,'(a)') 'beta_dfa.dat is opened!'
+
+ read(iodfa,'(a)') buffer
+ read(iodfa,*) itmp
+ read(iodfa,*) beta_inc
+
+ do i=1,itmp
+ read(iodfa,*) ica1, iitmp
+ do j=1,itmp
+ read(iodfa,*) wtmp
+ wshet(i,j) = wtmp
+c write(*,*) 'BETA:',i,j,wtmp,wshet(i,j)
+ enddo
+ enddo
+
+ close(iodfa)
+C END OF BETA RESTRAINT
+
+ return
+ END
+
+ subroutine edfad(edfadis)
+
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.DFA'
+
+ double precision edfadis
+ integer i, iatm1, iatm2,idiff
+ double precision ckk, sckk,dist,texp
+ double precision jix,jiy,jiz,ep,fp,scc
+
+ edfadis=0
+ gdfad=0.0d0
+
+ do i=1, idfadis
+
+ iatm1=idislis(1,i)+ishiftca
+ iatm2=idislis(2,i)+ishiftca
+ idiff = abs(iatm1-iatm2)
+
+ JIX=c(1,iatm2)-c(1,iatm1)
+ JIY=c(2,iatm2)-c(2,iatm1)
+ JIZ=c(3,iatm2)-c(3,iatm1)
+ DIST=SQRT(JIX*JIX+JIY*JIY+JIZ*JIZ)
+
+ ckk=ck(idiff)
+ sckk=sck(idiff)
+
+ scc = 0.0d0
+ ep = 0.0d0
+ fp = 0.0d0
+
+ do j=1,idisnum(i)
+
+ dd = dist-fdist(i,j)
+ dtmp = dd*dd/ckk
+ if (dtmp.ge.15.0d0) then
+ texp = 0.0d0
+ else
+c texp = dfaexp( idint(dtmp*1000)+1 )/sckk
+ texp = exp(-dtmp)/sckk
+ endif
+
+ ep=ep+sccdist(i,j)*texp
+ fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
+ scc=scc+sccdist(i,j)
+C write(*,'(2i8,6f12.5)') i, j, dist,
+C & fdist(i,j), ep, fp, sccdist(i,j), scc
+
+ enddo
+
+ ep = -ep/scc
+ fp = fp/scc
+
+
+c IF(ABS(EP).lt.1.0d-20)THEN
+c EP=0.0D0
+c ENDIF
+c IF (ABS(FP).lt.1.0d-20) THEN
+c FP=0.0D0
+c ENDIF
+
+ edfadis=edfadis+ep*dis_inc*wwdist
+
+ gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
+ gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
+ gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
+
+ gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
+ gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
+ gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
+
+ enddo
+
+ return
+ end
+
+ subroutine edfat(edfator)
+C DFA torsion angle
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.DFA'
+
+ integer i,j,ii,iii
+ integer iatom(5)
+ double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
+ double precision cwidth, cwidth2
+ PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
+
+ edfator= 0.0d0
+ enephi = 0.0d0
+ enethe = 0.0d0
+ gdfat(:,:) = 0.0d0
+
+C START OF PHI ANGLE
+ do i=1, idfaphi
+
+ aphi = 0.0d0
+ do iii=1,5
+ iatom(iii)=iphilis(iii,i)+ishiftca
+ enddo
+
+C ANGLE VECTOR CALCULTION
+ RIX=C(1,IATOM(2))-C(1,IATOM(1))
+ RIY=C(2,IATOM(2))-C(2,IATOM(1))
+ RIZ=C(3,IATOM(2))-C(3,IATOM(1))
+
+ RIPX=C(1,IATOM(3))-C(1,IATOM(2))
+ RIPY=C(2,IATOM(3))-C(2,IATOM(2))
+ RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
+
+ RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
+ RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
+ RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
+
+ RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
+ RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
+ RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
+
+ GIX=RIY*RIPZ-RIZ*RIPY
+ GIY=RIZ*RIPX-RIX*RIPZ
+ GIZ=RIX*RIPY-RIY*RIPX
+
+ GIPX=RIPY*RIPPZ-RIPZ*RIPPY
+ GIPY=RIPZ*RIPPX-RIPX*RIPPZ
+ GIPZ=RIPX*RIPPY-RIPY*RIPPX
+
+ CIPX=C(1,IATOM(3))-C(1,IATOM(1))
+ CIPY=C(2,IATOM(3))-C(2,IATOM(1))
+ CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
+
+ CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
+ CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
+ CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
+
+ CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
+ CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
+ CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
+
+ DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
+ DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
+ DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
+ DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
+
+C END OF ANGLE VECTOR CALCULTION
+
+ TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
+ APHI(1)=TDOT/(DGI*DRIPP)
+ TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
+ APHI(2)=TDOT/(DGIP*DRIP3)
+
+ ephi = 0.0d0
+ tfphi1=0.0d0
+ tfphi2=0.0d0
+ scc=0.0d0
+
+ do j=1, iphinum(i)
+ DDPS1=APHI(1)-FPHI1(i,j)
+ DDPS2=APHI(2)-FPHI2(i,j)
+
+ DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
+
+ if (dtmp.ge.15.0d0) then
+ ps_tmp = 0.0d0
+ else
+c ps_tmp = dfaexp(idint(dtmp*1000)+1)
+ ps_tmp = exp(-dtmp)
+ endif
+
+ ephi=ephi+sccphi(i,j)*ps_tmp
+
+ tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
+ tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
+
+ scc=scc+sccphi(i,j)
+C write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
+C & aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
+ ENDDO
+
+ ephi=-ephi/scc*phi_inc*wwangle
+ tfphi1=tfphi1/scc*phi_inc*wwangle
+ tfphi2=tfphi2/scc*phi_inc*wwangle
+
+ IF (ABS(EPHI).LT.1d-20) THEN
+ EPHI=0.0D0
+ ENDIF
+ IF (ABS(TFPHI1).LT.1d-20) THEN
+ TFPHI1=0.0D0
+ ENDIF
+ IF (ABS(TFPHI2).LT.1d-20) THEN
+ TFPHI2=0.0D0
+ ENDIF
+
+C FORCE DIRECTION CALCULATION
+ TDX(1:5)=0.0D0
+ TDY(1:5)=0.0D0
+ TDZ(1:5)=0.0D0
+
+ DM1=1.0d0/(DGI*DRIPP)
+
+ GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
+ DM2=GIRPP/(DGI**3*DRIPP)
+ DM3=GIRPP/(DGI*DRIPP**3)
+
+ DM4=1.0d0/(DGIP*DRIP3)
+
+ GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
+ DM5=GIRP3/(DGIP**3*DRIP3)
+ DM6=GIRP3/(DGIP*DRIP3**3)
+C FIRST ATOM BY PHI1
+ TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
+ & +( GIZ* RIPY- GIY* RIPZ)*DM2
+ TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
+ & +( GIX* RIPZ- GIZ* RIPX)*DM2
+ TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
+ & +( GIY* RIPX- GIX* RIPY)*DM2
+ TDX(1)=TDX(1)*TFPHI1
+ TDY(1)=TDY(1)*TFPHI1
+ TDZ(1)=TDZ(1)*TFPHI1
+C SECOND ATOM BY PHI1
+ TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
+ & -(CIPY*GIZ-CIPZ*GIY)*DM2
+ TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
+ & -(CIPZ*GIX-CIPX*GIZ)*DM2
+ TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
+ & -(CIPX*GIY-CIPY*GIX)*DM2
+ TDX(2)=TDX(2)*TFPHI1
+ TDY(2)=TDY(2)*TFPHI1
+ TDZ(2)=TDZ(2)*TFPHI1
+C SECOND ATOM BY PHI2
+ TDX(2)=TDX(2)+
+ & ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
+ & +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
+ TDY(2)=TDY(2)+
+ & ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
+ & +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
+ TDZ(2)=TDZ(2)+
+ & ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
+ & +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
+C THIRD ATOM BY PHI1
+ TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
+ & -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
+ TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
+ & -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
+ TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
+ & -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
+ TDX(3)=TDX(3)*TFPHI1
+ TDY(3)=TDY(3)*TFPHI1
+ TDZ(3)=TDZ(3)*TFPHI1
+C THIRD ATOM BY PHI2
+ TDX(3)=TDX(3)+
+ & ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
+ & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
+ TDY(3)=TDY(3)+
+ & ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
+ & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
+ TDZ(3)=TDZ(3)+
+ & ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
+ & -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
+C FOURTH ATOM BY PHI1
+ TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
+ TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
+ TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
+C FOURTH ATOM BY PHI2
+ TDX(4)=TDX(4)+
+ & ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
+ & -( GIPY*RIPZ-RIPY*GIPZ)*DM5
+ & + RIP3X*DM6)*TFPHI2
+ TDY(4)=TDY(4)+
+ & ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
+ & -( GIPZ*RIPX-RIPZ*GIPX)*DM5
+ & + RIP3Y*DM6)*TFPHI2
+ TDZ(4)=TDZ(4)+
+ & ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
+ & -( GIPX*RIPY-RIPX*GIPY)*DM5
+ & + RIP3Z*DM6)*TFPHI2
+C FIFTH ATOM BY PHI2
+ TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
+ TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
+ TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
+C END OF FORCE DIRECTION
+c force calcuation
+ DO II=1,5
+ gdfat(1,IATOM(II))=gdfat(1,IATOM(II))+TDX(II)
+ gdfat(2,IATOM(II))=gdfat(2,IATOM(II))+TDY(II)
+ gdfat(3,IATOM(II))=gdfat(3,IATOM(II))+TDZ(II)
+ ENDDO
+c energy calculation
+ enephi = enephi + ephi
+c end of single assignment statement
+ ENDDO
+C END OF PHI RESTRAINT
+
+C START OF THETA ANGLE
+ do i=1, idfathe
+
+ athe = 0.0d0
+ do iii=1,5
+ iatom(iii)=ithelis(iii,i)+ishiftca
+ enddo
+
+
+C ANGLE VECTOR CALCULTION
+ RIX=C(1,IATOM(2))-C(1,IATOM(1))
+ RIY=C(2,IATOM(2))-C(2,IATOM(1))
+ RIZ=C(3,IATOM(2))-C(3,IATOM(1))
+
+ RIPX=C(1,IATOM(3))-C(1,IATOM(2))
+ RIPY=C(2,IATOM(3))-C(2,IATOM(2))
+ RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
+
+ RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
+ RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
+ RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
+
+ RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
+ RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
+ RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
+
+ GIX=RIY*RIPZ-RIZ*RIPY
+ GIY=RIZ*RIPX-RIX*RIPZ
+ GIZ=RIX*RIPY-RIY*RIPX
+
+ GIPX=RIPY*RIPPZ-RIPZ*RIPPY
+ GIPY=RIPZ*RIPPX-RIPX*RIPPZ
+ GIPZ=RIPX*RIPPY-RIPY*RIPPX
+
+ GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
+ GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
+ GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
+
+ CIPX=C(1,IATOM(3))-C(1,IATOM(1))
+ CIPY=C(2,IATOM(3))-C(2,IATOM(1))
+ CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
+
+ CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
+ CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
+ CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
+
+ CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
+ CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
+ CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
+
+ DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
+ DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
+ DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
+ DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
+ DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
+C END OF ANGLE VECTOR CALCULTION
+
+ TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
+ ATHE(1)=TDOT/(DGI*DGIP)
+ TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
+ ATHE(2)=TDOT/(DGIP*DGIPP)
+
+ ETHE=0.0D0
+ TFTHE1=0.0D0
+ TFTHE2=0.0D0
+ SCC=0.0D0
+ TH_TMP=0.0d0
+
+ do j=1,ithenum(i)
+ ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
+ ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
+ dtmp= (ddth1**2+ddth2**2)/cwidth2
+ if ( dtmp .ge. 15.0d0) then
+ th_tmp = 0.0d0
+ else
+c th_tmp = dfaexp ( idint(dtmp*1000)+1 )
+ th_tmp = exp(-dtmp)
+ end if
+
+ ethe=ethe+sccthe(i,j)*th_tmp
+
+ tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
+ tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
+ scc=scc+sccthe(i,j)
+C write(*,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
+C & athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
+ enddo
+
+ ethe=-ethe/scc*the_inc*wwangle
+ tfthe1=tfthe1/scc*the_inc*wwangle
+ tfthe2=tfthe2/scc*the_inc*wwangle
+
+ IF (ABS(ETHE).LT.TENM20) THEN
+ ETHE=0.0D0
+ ENDIF
+ IF (ABS(TFTHE1).LT.TENM20) THEN
+ TFTHE1=0.0D0
+ ENDIF
+ IF (ABS(TFTHE2).LT.TENM20) THEN
+ TFTHE2=0.0D0
+ ENDIF
+
+ TDX(1:5)=0.0D0
+ TDY(1:5)=0.0D0
+ TDZ(1:5)=0.0D0
+
+ DM1=1.0d0/(DGI*DGIP)
+ DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
+ DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
+
+ DM4=1.0d0/(DGIP*DGIPP)
+ DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
+ DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
+
+C FIRST ATOM BY THETA1
+ TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
+ & -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
+ TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
+ & -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
+ TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
+ & -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
+C SECOND ATOM BY THETA1
+ TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
+ & -(CIPY*GIZ-CIPZ*GIY)*DM2
+ & +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
+ TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
+ & -(CIPZ*GIX-CIPX*GIZ)*DM2
+ & +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
+ TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
+ & -(CIPX*GIY-CIPY*GIX)*DM2
+ & +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
+C SECOND ATOM BY THETA2
+ TDX(2)=TDX(2)+
+ & ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
+ & -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
+ TDY(2)=TDY(2)+
+ & ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
+ & -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
+ TDZ(2)=TDZ(2)+
+ & ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
+ & -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
+C THIRD ATOM BY THETA1
+ TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
+ & -(GIY*RIZ-GIZ*RIY)*DM2
+ & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
+ TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
+ & -(GIZ*RIX-GIX*RIZ)*DM2
+ & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
+ TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
+ & -(GIX*RIY-GIY*RIX)*DM2
+ & -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
+C THIRD ATOM BY THETA2
+ TDX(3)=TDX(3)+
+ & ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
+ & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
+ & +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
+ TDY(3)=TDY(3)+
+ & ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
+ & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
+ & +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
+ TDZ(3)=TDZ(3)+
+ & ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
+ & -(CIPPX*GIPY-CIPPY*GIPX)*DM5
+ & +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
+C FOURTH ATOM BY THETA1
+ TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
+ & -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
+ TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
+ & -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
+ TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
+ & -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
+C FOURTH ATOM BY THETA2
+ TDX(4)=TDX(4)+
+ & ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
+ & -(GIPY*RIPZ-GIPZ*RIPY)*DM5
+ & -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
+ TDY(4)=TDY(4)+
+ & ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
+ & -(GIPZ*RIPX-GIPX*RIPZ)*DM5
+ & -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
+ TDZ(4)=TDZ(4)+
+ & ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
+ & -(GIPX*RIPY-GIPY*RIPX)*DM5
+ & -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
+C FIFTH ATOM BY THETA2
+ TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
+ & -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
+ TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
+ & -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
+ TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
+ & -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
+C !! END OF FORCE DIRECTION!!!!
+ DO II=1,5
+ gdfat(1,iatom(II))=gdfat(1,iatom(II))+TDX(II)
+ gdfat(2,iatom(II))=gdfat(2,iatom(II))+TDY(II)
+ gdfat(3,iatom(II))=gdfat(3,iatom(II))+TDZ(II)
+ ENDDO
+C energy calculation
+ enethe = enethe + ethe
+ ENDDO
+
+ edfator = enephi + enethe
+
+ RETURN
+ END
+
+ subroutine edfan(edfanei)
+C DFA neighboring CA restraint
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.DFA'
+
+ integer i,j,imin
+ integer kshnum, n1atom
+
+ double precision enenei,tmp_n
+ double precision pai,hpai
+ double precision jix,jiy,jiz,ndiff,snorm_nei
+ double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
+ double precision dr,dr2,half,ntmp,dtmp
+
+ parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
+ parameter(pai=3.14159265358979323846D0)
+ parameter(hpai=1.5707963267948966D0)
+ parameter(snorm_nei=0.886226925452758D0)
+
+ edfanei = 0.0d0
+ enenei = 0.0d0
+ gdfan = 0.0d0
+
+c print*, 's1:', s1(:)
+c print*, 's2:', s2(:)
+
+ do i=1, idfanei
+
+ kshnum=kshell(i)
+ n1atom=ineilis(i)+ishiftca
+C write(*,*) 'kshnum,n1atom:', kshnum, n1atom
+
+ tmp_n=0.0d0
+ ftmp=0.0d0
+ dnei=0.0d0
+ dist=0.0d0
+ t1dx=0.0d0
+ t1dy=0.0d0
+ t1dz=0.0d0
+ t2dx=0.0d0
+ t2dy=0.0d0
+ t2dz=0.0d0
+
+ do j = ishiftca+1, ilastca
+
+ if (n1atom.eq.j) cycle
+
+ jix=c(1,j)-c(1,n1atom)
+ jiy=c(2,j)-c(2,n1atom)
+ jiz=c(3,j)-c(3,n1atom)
+ dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
+
+c write(*,*) n1atom, j, dist
+
+ if(kshnum.ne.1)then
+ if (dist.lt.s1(kshnum).and.
+ & dist.gt.s2(kshnum-1)) then
+
+ tmp_n=tmp_n+1.0d0
+
+c write(*,*) 'case1:',tmp_n
+
+ t1dx=t1dx+0.0d0
+ t1dy=t1dy+0.0d0
+ t1dz=t1dz+0.0d0
+ t2dx(j)=0.0d0
+ t2dy(j)=0.0d0
+ t2dz(j)=0.0d0
+
+ elseif(dist.ge.s1(kshnum).and.
+ & dist.le.s2(kshnum)) then
+
+ dnei=(dist-s1(kshnum))/dr2*pai
+ tmp_n=tmp_n + half*(1+cos(dnei))
+c write(*,*) 'case2:',tmp_n
+ ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
+c center atom
+ t1dx=t1dx+jix*ftmp
+ t1dy=t1dy+jiy*ftmp
+ t1dz=t1dz+jiz*ftmp
+c neighbor atoms
+ t2dx(j)=-jix*ftmp
+ t2dy(j)=-jiy*ftmp
+ t2dz(j)=-jiz*ftmp
+c
+ elseif(dist.ge.s1(kshnum-1).and.
+ & dist.le.s2(kshnum-1)) then
+ dnei=(dist-s1(kshnum-1))/dr2*pai
+ tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
+c write(*,*) 'case3:',tmp_n
+ ftmp = hpai*sin(dnei)/dr2/dist
+c center atom
+ t1dx=t1dx+jix*ftmp
+ t1dy=t1dy+jiy*ftmp
+ t1dz=t1dz+jiz*ftmp
+c neighbor atoms
+ t2dx(j)=-jix*ftmp
+ t2dy(j)=-jiy*ftmp
+ t2dz(j)=-jiz*ftmp
+
+ endif
+
+ elseif(kshnum.eq.1) then
+
+ if(dist.lt.s1(kshnum))then
+
+ tmp_n=tmp_n+1.0d0
+c write(*,*) 'case4:',tmp_n
+ t1dx=t1dx+0.0d0
+ t1dy=t1dy+0.0d0
+ t1dz=t1dz+0.0d0
+ t2dx(j)=0.0d0
+ t2dy(j)=0.0d0
+ t2dz(j)=0.0d0
+
+ elseif(dist.ge.s1(kshnum).and.
+ & dist.le.s2(kshnum))then
+
+ dnei=(dist-s1(kshnum))/dr2*pai
+ tmp_n=tmp_n + half*(1+cos(dnei))
+c write(*,*) 'case5:',tmp_n
+ ftmp = -hpai*sin(dnei)/dr2/dist
+c center atom
+ t1dx=t1dx+jix*ftmp
+ t1dy=t1dy+jiy*ftmp
+ t1dz=t1dz+jiz*ftmp
+c neighbor atoms
+ t2dx(j)=-jix*ftmp
+ t2dy(j)=-jiy*ftmp
+ t2dz(j)=-jiz*ftmp
+
+ endif
+ endif
+ enddo
+
+ scc=0.0d0
+ enei=0.0d0
+ tmp_fnei=0.0d0
+ ndiff=0.0d0
+
+ do imin=1,ineinum(i)
+
+ ndiff = tmp_n-fnei(i,imin)
+ dtmp = ndiff*ndiff
+
+ if (dtmp.ge.15.0d0) then
+ ntmp = 0.0d0
+ else
+c ntmp = dfaexp( idint(dtmp*1000) + 1 )
+ ntmp = exp(-dtmp)
+ end if
+
+ enei=enei+sccnei(i,imin)*ntmp
+ tmp_fnei=tmp_fnei-
+ & sccnei(i,imin)*ntmp*ndiff*2.0d0
+ scc=scc+sccnei(i,imin)
+
+c write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
+c & fnei(i,imin),sccnei(i,imin),enei,scc
+ enddo
+
+ enei=-enei/scc*snorm_nei*nei_inc*wwnei
+ tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
+
+c if (abs(enei).lt.1.0d-20)then
+c enei=0.0d0
+c endif
+c if (abs(tmp_fnei).lt.1.0d-20) then
+c tmp_fnei=0.0d0
+c endif
+
+c force calculation
+ t1dx=t1dx*tmp_fnei
+ t1dy=t1dy*tmp_fnei
+ t1dz=t1dz*tmp_fnei
+
+ do j=ishiftca+1,ilastca
+ t2dx(j)=t2dx(j)*tmp_fnei
+ t2dy(j)=t2dy(j)*tmp_fnei
+ t2dz(j)=t2dz(j)*tmp_fnei
+ enddo
+
+ gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
+ gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
+ gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
+
+ do j=ishiftca+1,ilastca
+ gdfan(1,j)=gdfan(1,j)+t2dx(j)
+ gdfan(2,j)=gdfan(2,j)+t2dy(j)
+ gdfan(3,j)=gdfan(3,j)+t2dz(j)
+ enddo
+c energy calculation
+
+ enenei=enenei+enei
+
+ enddo
+
+ edfanei=enenei
+
+ return
+ end
+
+ subroutine edfab(edfabeta)
+
+ implicit real*8 (a-h,o-z)
+
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.DFA'
+
+ real*8 PAI
+ parameter(PAI=3.14159265358979323846D0)
+ parameter (maxca=800)
+C sheet variables
+ real*8 bx(maxres),by(maxres),bz(maxres)
+ real*8 vbet(maxres,maxres)
+ real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
+ real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
+ real*8 vbeta,vbetp,vbetm
+ real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ & c00,s00,ulnex,dnex
+ real*8 dp45,dm45,w_beta
+
+ real*8 cph(maxca),cth(maxca)
+ real*8 atx(maxca),aty(maxca),atz(maxca)
+ real*8 atmx(maxca),atmy(maxca),atmz(maxca)
+ real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
+ real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
+ real*8 sth(maxca)
+ real*8 astx(maxca),asty(maxca),astz(maxca)
+ real*8 astmx(maxca),astmy(maxca),astmz(maxca)
+ real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
+ real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
+
+ real*8 atxnum(maxca),atynum(maxca),atznum(maxca),
+ & astxnum(maxca),astynum(maxca),astznum(maxca),
+ & atmxnum(maxca),atmynum(maxca),atmznum(maxca),
+ & astmxnum(maxca),astmynum(maxca),astmznum(maxca),
+ & atmmxnum(maxca),atmmynum(maxca),atmmznum(maxca),
+ & astmmxnum(maxca),astmmynum(maxca),astmmznum(maxca),
+ & atm3xnum(maxca),atm3ynum(maxca),atm3znum(maxca),
+ & astm3xnum(maxca),astm3ynum(maxca),astm3znum(maxca),
+ & cth_orig(maxca),sth_orig(maxca)
+
+ common /sheca/ bx,by,bz
+ common /shee/ vbeta,vbet,vbetp,vbetm
+ common /shetf/ shetfx,shetfy,shetfz
+ common /shef/ shefx, shefy, shefz
+ common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ & c00,s00,ulnex,dnex
+ common /sheconst/ dp45,dm45,w_beta
+
+ common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
+ $ atmmz,atm3x,atm3y,atm3z
+ common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
+ $ astmmz,astm3x,astm3y,astm3z
+
+ common /coscos/ cph,cth
+ common /sinsin/ sth
+
+C End of sheet variables
+
+ integer i,j
+ double precision enebet
+
+ enebet=0.0d0
+ bx=0.0d0;by=0.0d0;bz=0.0d0
+ shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
+
+ gdfab=0.0d0
+
+ do i=ishiftca+1,ilastca
+ bx(i-ishiftca)=c(1,i)
+ by(i-ishiftca)=c(2,i)
+ bz(i-ishiftca)=c(3,i)
+ enddo
+
+c do i=1,ilastca-ishiftca
+c read(99,*) bx(i),by(i),bz(i)
+c enddo
+c close(99)
+
+ dca=0.25d0**2
+ dshe=0.3d0**2
+ ULHB=5.0D0
+ ULDHB=5.0D0
+ ULNEX=COS(60.0D0/180.0D0*PAI)
+
+ DLHB=1.0D0
+ DLDHB=1.0D0
+
+ DNEX=0.3D0**2
+
+ C00=COS((1.0D0+10.0D0/180.0D0)*PAI)
+ S00=SIN((1.0D0+10.0D0/180.0D0)*PAI)
+
+ W_BETA=0.5D0
+ DP45=W_BETA
+ DM45=W_BETA
+
+C END OF INITIALIZATION
+
+ nca=ilastca-ishiftca
+
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+
+c end of sheet energy and force
+
+ do j=1,nca
+ shetfx(j)=shetfx(j)*beta_inc
+ shetfy(j)=shetfy(j)*beta_inc
+ shetfz(j)=shetfz(j)*beta_inc
+c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
+ enddo
+
+ vbeta=vbeta*beta_inc
+ enebet=vbeta
+ edfabeta=enebet
+
+ do j=1,nca
+ gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
+ gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
+ gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
+ enddo
+
+#ifdef DEBUG1
+ do j=1,nca
+ write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
+ enddo
+
+
+ gdfab=0
+ dinc=0.001
+ do j=1,nca
+ cth_orig(j)=cth(j)
+ sth_orig(j)=sth(j)
+ enddo
+
+ do j=1,nca
+
+ bx(j)=bx(j)+dinc
+ call angvectors(nca)
+ bx(j)=bx(j)-2*dinc
+ call angvectors(nca)
+ atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
+ astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
+ if (j.gt.1) then
+ atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
+ astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
+ endif
+ if (j.gt.2) then
+ atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
+ astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
+ endif
+ if (j.gt.3) then
+ atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
+ astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
+ endif
+ bx(j)=bx(j)+dinc
+ by(j)=by(j)+dinc
+ call angvectors(nca)
+ by(j)=by(j)-2*dinc
+ call angvectors(nca)
+ by(j)=by(j)+dinc
+ atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
+ astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
+ if (j.gt.1) then
+ atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
+ astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
+ endif
+ if (j.gt.2) then
+ atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
+ astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
+ endif
+ if (j.gt.3) then
+ atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
+ astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
+ endif
+
+ bz(j)=bz(j)+dinc
+ call angvectors(nca)
+ bz(j)=bz(j)-2*dinc
+ call angvectors(nca)
+ bz(j)=bz(j)+dinc
+
+ atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
+ astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
+ if (j.gt.1) then
+ atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
+ astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
+ endif
+ if (j.gt.2) then
+ atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
+ astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
+ endif
+ if (j.gt.3) then
+ atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
+ astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
+ endif
+
+ enddo
+
+ do i=1,nca
+ write (*,'(2i5,a2,6f10.5)')
+ & i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
+ & astxnum(i),astx(i),astxnum(i)/astx(i),
+ & i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
+ & astynum(i),asty(i),astynum(i)/asty(i),
+ & i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
+ & astznum(i),astz(i),astznum(i)/astz(i),
+ & i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
+ & astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
+ & i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
+ & astmynum(i),astmy(i),astmynum(i)/astmy(i),
+ & i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
+ & astmznum(i),astmz(i),astmznum(i)/astmz(i),
+ & i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
+ & astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
+ & i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
+ & astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
+ & i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
+ & astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
+ & i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
+ & astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
+ & i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
+ & astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
+ & i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
+ & astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
+ & i,0," ",cth_orig(i),sth_orig(i)
+ enddo
+
+
+ gdfab=0
+ dinc=0.001
+
+ do j=1,nca
+
+ bx(j)=bx(j)+dinc
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+ vbeta1=vbeta*beta_inc
+ bx(j)=bx(j)-2*dinc
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+ vbeta2=vbeta*beta_inc
+ gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
+ bx(j)=bx(j)+dinc
+
+ by(j)=by(j)+dinc
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+ vbeta1=vbeta*beta_inc
+ by(j)=by(j)-2*dinc
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+ vbeta2=vbeta*beta_inc
+ gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
+ by(j)=by(j)+dinc
+
+ bz(j)=bz(j)+dinc
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+ vbeta1=vbeta*beta_inc
+ bz(j)=bz(j)-2*dinc
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+ vbeta2=vbeta*beta_inc
+ gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
+ bz(j)=bz(j)+dinc
+
+
+ enddo
+
+
+ call angvectors(nca)
+ call sheetforce(nca,wshet)
+ do j=1,nca
+ shetfx(j)=shetfx(j)*beta_inc
+ shetfy(j)=shetfy(j)*beta_inc
+ shetfz(j)=shetfz(j)*beta_inc
+ enddo
+
+
+ write(*,*) 'xyz analytical and numerical gradient'
+ do j=1,nca
+ write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
+ & ,(-gdfab(i,j),i=1,3)
+ enddo
+
+ do j=1,nca
+ write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
+ & shetfy(j)/gdfab(2,j),
+ & shetfz(j)/gdfab(3,j)
+ enddo
+
+ stop
+#endif
+
+ return
+ end
+C-------------------------------------------------------------------------------
+ subroutine angvectors(nca)
+c implicit real*4(a-h,o-z)
+ implicit none
+ integer nca
+ integer maxca
+ parameter(maxca=800)
+ real*8 pai,zero
+ parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
+
+ real*8 bx(maxca),by(maxca),bz(maxca)
+ real*8 dis(maxca,maxca)
+ real*8 apx(maxca),apy(maxca),apz(maxca)
+ real*8 apmx(maxca),apmy(maxca),apmz(maxca)
+ real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
+ real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
+ real*8 atx(maxca),aty(maxca),atz(maxca)
+ real*8 atmx(maxca),atmy(maxca),atmz(maxca)
+ real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
+ real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
+ real*8 astx(maxca),asty(maxca),astz(maxca)
+ real*8 astmx(maxca),astmy(maxca),astmz(maxca)
+ real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
+ real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
+ real*8 sth(maxca)
+ real*8 cph(maxca),cth(maxca)
+ real*8 ulcos(maxca)
+ real*8 p,c
+ integer i, ip, ipp, ip3, j
+ real*8 rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
+ real*8 rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
+ real*8 gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
+ real*8 cix, ciy, ciz, cipx, cipy, cipz
+ real*8 gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
+ real*8 d10, d11, d12, d13, d20, d21, d22, d23, d24
+ real*8 d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
+ real*8 d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
+ real*8 dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
+ real*8 dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
+ real*8 g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
+ real*8 gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
+ real*8 gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
+ real*8 gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
+ real*8 gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
+ real*8 grpp,gx,gy,gz
+ real*8 rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
+ real*8 sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
+ integer inb,nmax,iselect
+
+ common /sheca/ bx,by,bz
+ common /difvec/ rx, ry, rz
+ common /ulang/ ulcos
+ common /phys1/ inb,nmax,iselect
+ common /phys4/ p,c
+ common /kyori2/ dis
+ common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
+ & apmmz,apm3x,apm3y,apm3z
+ common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
+ & atmmz,atm3x,atm3y,atm3z
+ common /coscos/ cph,cth
+ common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
+ & astmmz,astm3x,astm3y,astm3z
+ common /sinsin/ sth
+C-------------------------------------------------------------------------------
+c write(*,*) 'inside angvectors'
+C initialize
+ p=0.1d0
+ c=1.0d0
+ inb=nca
+ cph=zero; cth=zero; sth=zero
+ apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
+ apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
+ atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
+ atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
+ astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
+ astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
+ astm3z=zero
+C end of initialize
+C r[x,y,z] calc and distance calculation
+ rx=zero;ry=zero;rz=zero
+
+ do i=1,inb
+ do j=1,inb
+ rx(i,j)=bx(j)-bx(i)
+ ry(i,j)=by(j)-by(i)
+ rz(i,j)=bz(j)-bz(i)
+ dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
+c write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
+c write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
+c write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
+c write(*,*) 'dis(i,j):',i,j,dis(i,j)
+ enddo
+ enddo
+c end of r[x,y,z] calc
+C cos calc
+ do i=1,inb-2
+ ip=i+1
+ ipp=i+2
+
+ if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
+ ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
+ $ +rz(i,ip)*rz(ip,ipp)
+ ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
+ endif
+ enddo
+c end of virtual bond angle
+c write(*,*) 'inside angvectors1'
+crc do i=1,inb-3
+ do i=1,inb
+ ip=i+1
+ ipp=i+2
+ ip3=i+3
+ rix=bx(ip)-bx(i)
+ riy=by(ip)-by(i)
+ riz=bz(ip)-bz(i)
+ ripx=bx(ipp)-bx(ip)
+ ripy=by(ipp)-by(ip)
+ ripz=bz(ipp)-bz(ip)
+ rippx=bx(ip3)-bx(ipp)
+ rippy=by(ip3)-by(ipp)
+ rippz=bz(ip3)-bz(ipp)
+
+ gx=riy*ripz-riz*ripy
+ gy=riz*ripx-rix*ripz
+ gz=rix*ripy-riy*ripx
+ gpx=ripy*rippz-ripz*rippy
+ gpy=ripz*rippx-ripx*rippz
+ gpz=ripx*rippy-ripy*rippx
+ gpcrp_x=gpy*ripz-gpz*ripy
+ gpcrp_y=gpz*ripx-gpx*ripz
+ gpcrp_z=gpx*ripy-gpy*ripx
+ d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
+ gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
+ & -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
+
+ if(i.ge.2) then
+ rimx=bx(i)-bx(i-1)
+ rimy=by(i)-by(i-1)
+ rimz=bz(i)-bz(i-1)
+ gmx=rimy*riz-rimz*riy
+ gmy=rimz*rix-rimx*riz
+ gmz=rimx*riy-rimy*rix
+ dgm=sqrt(gmx**2+gmy**2+gmz**2)
+ dgm3=dgm**3
+ ggm=gmx*gx+gmy*gy+gmz*gz
+ gmrp=gmx*ripx+gmy*ripy+gmz*ripz
+ drim=dis(i-1,i)
+ drim3=drim**3
+ gcr_x=gy*riz-gz*riy
+ gcr_y=gz*rix-gx*riz
+ gcr_z=gx*riy-gy*rix
+ d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
+ d_gcr3=d_gcr**3
+ gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
+ & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
+ endif
+c write(*,*) 'inside angvectors2'
+ if(i.ge.3) then
+ rimmx=bx(i-1)-bx(i-2)
+ rimmy=by(i-1)-by(i-2)
+ rimmz=bz(i-1)-bz(i-2)
+ drimm=dis(i-2,i-1)
+ gmmx=rimmy*rimz-rimmz*rimy
+ gmmy=rimmz*rimx-rimmx*rimz
+ gmmz=rimmx*rimy-rimmy*rimx
+ dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
+ dgmm3=dgmm**3
+ gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
+ gmmr=gmmx*rix+gmmy*riy+gmmz*riz
+ gmcrim_x=gmy*rimz-gmz*rimy
+ gmcrim_y=gmz*rimx-gmx*rimz
+ gmcrim_z=gmx*rimy-gmy*rimx
+ d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
+ d_gmcrim3=d_gmcrim**3
+ gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
+ & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
+ endif
+
+ if(i.ge.4) then
+ rim3x=bx(i-2)-bx(i-3)
+ rim3y=by(i-2)-by(i-3)
+ rim3z=bz(i-2)-bz(i-3)
+ g3x=rim3y*rimmz-rim3z*rimmy
+ g3y=rim3z*rimmx-rim3x*rimmz
+ g3z=rim3x*rimmy-rim3y*rimmx
+ dg30=sqrt(g3x**2+g3y**2+g3z**2)
+ g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
+ g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
+cc**********************************************************************
+ gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
+ gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
+ gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
+ d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
+ d_gmmcrimm3=d_gmmcrimm**3
+ gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
+ & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
+ endif
+
+ dri=dis(i,i+1)
+ drip=dis(i+1,i+2)
+ dripp=dis(i+2,i+3)
+ dri3=dri**3
+ dg=sqrt(gx**2+gy**2+gz**2)
+ dgp=sqrt(gpx**2+gpy**2+gpz**2)
+ dg3=dg**3
+
+ ggp=gx*gpx+gy*gpy+gz*gpz
+ grpp=gx*rippx+gy*rippy+gz*rippz
+
+ if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
+ & .and.d_gpcrp.gt.0.0D0) then
+ cph(i)=grpp/dg/dripp
+ cth(i)=ggp/dg/dgp
+ sth(i)=gpcrp__g/d_gpcrp/dg
+ else
+c
+ cph(i)=1.0D0
+ cth(i)=1.0D0
+ sth(i)=0.0D0
+ endif
+
+c write(*,*) 'inside angvectors3'
+
+ if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
+ & .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
+ d10=1.0D0/(dg*dgp)
+ d11=ggp/(dg3*dgp)
+ d12=1.0D0/(dg*dripp)
+ d13=grpp/(dg3*dripp)
+ sd10=1.0D0/(d_gpcrp*dg)
+ sd11=gpcrp__g/(d_gpcrp*dg3)
+ else
+ d10=0.0D0
+ d11=0.0D0
+ d12=0.0D0
+ d13=0.0D0
+ sd10=0.0D0
+ sd11=0.0D0
+ endif
+
+ atx(i)=(ripz*gpy-ripy*gpz)*d10
+ & -(gy*ripz-gz*ripy)*d11
+ aty(i)=(ripx*gpz-ripz*gpx)*d10
+ & -(gz*ripx-gx*ripz)*d11
+ atz(i)=(ripy*gpx-ripx*gpy)*d10
+ & -(gx*ripy-gy*ripx)*d11
+ astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
+ & +ripy*gpy*ripx-gpx*ripz**2)
+ & -sd11*(gy*ripz-gz*ripy)
+ asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
+ & -gpy*ripx**2+gpz*ripy*ripz)
+ & -sd11*(-gx*ripz+gz*ripx)
+ astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
+ & -gpz*ripy**2+ripz*gpx*ripx)
+ & -sd11*(gx*ripy-gy*ripx)
+ apx(i)=(ripz*rippy-ripy*rippz)*d12
+ & -(gy*ripz-gz*ripy)*d13
+ apy(i)=(ripx*rippz-ripz*rippx)*d12
+ & -(gz*ripx-gx*ripz)*d13
+ apz(i)=(ripy*rippx-ripx*rippy)*d12
+ & -(gx*ripy-gy*ripx)*d13
+
+ if(i.ge.2) then
+ cix=bx(ip)-bx(i-1)
+ ciy=by(ip)-by(i-1)
+ ciz=bz(ip)-bz(i-1)
+ cipx=bx(ipp)-bx(i)
+ cipy=by(ipp)-by(i)
+ cipz=bz(ipp)-bz(i)
+ ripx=bx(ipp)-bx(ip)
+ ripy=by(ipp)-by(ip)
+ ripz=bz(ipp)-bz(ip)
+ if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
+ & .and.d_gcr3.gt.0.0D0) then
+ d20=1.0D0/(dg*dgm)
+ d21=ggm/(dgm3*dg)
+ d22=ggm/(dgm*dg3)
+ d23=1.0D0/(dgm*drip)
+ d24=gmrp/(dgm3*drip)
+ sd20=1.0D0/(d_gcr*dgm)
+ sd21=gcr__gm/(d_gcr3*dgm)
+ sd22=gcr__gm/(d_gcr*dgm3)
+ else
+ d20=0.0D0
+ d21=0.0D0
+ d22=0.0D0
+ d23=0.0D0
+ d24=0.0D0
+ sd20=0.0D0
+ sd21=0.0D0
+ sd22=0.0D0
+ endif
+ atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
+ & -(ciy*gmz-ciz*gmy)*d21
+ & +(ripy*gz-ripz*gy)*d22
+ atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
+ & -(ciz*gmx-cix*gmz)*d21
+ & +(ripz*gx-ripx*gz)*d22
+ atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
+ & -(cix*gmy-ciy*gmx)*d21
+ & +(ripx*gy-ripy*gx)*d22
+cc**********************************************************************
+ astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
+ & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
+ & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
+ & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
+ & +gcr_z*(-ripz*rix+gy))
+ & -sd22*(-gmy*ciz+gmz*ciy)
+
+ astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
+ & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
+ & +riz*ripz*gmy)
+ & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
+ & -gcr_z*(ripz*riy+gx))
+ & -sd22*(gmx*ciz-gmz*cix)
+
+ astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
+ & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
+ & -riz*gx*cix)
+ & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
+ & +gcr_z*(ripy*riy+ripx*rix))
+ & -sd22*(-gmx*ciy+gmy*cix)
+cc**********************************************************************
+ apmx(i)=(ciy*ripz-ripy*ciz)*d23
+ & -(ciy*gmz-ciz*gmy)*d24
+ apmy(i)=(ciz*ripx-ripz*cix)*d23
+ & -(ciz*gmx-cix*gmz)*d24
+ apmz(i)=(cix*ripy-ripx*ciy)*d23
+ & -(cix*gmy-ciy*gmx)*d24
+ endif
+
+ if(i.ge.3) then
+ if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
+ & .and.d_gmcrim3.gt.0.0D0) then
+ d30=1.0D0/(dgm*dgmm)
+ d31=gmmgm/(dgm3*dgmm)
+ d32=gmmgm/(dgm*dgmm3)
+ d33=1.0D0/(dgmm*dri)
+ d34=gmmr/(dgmm3*dri)
+ d35=gmmr/(dgmm*dri3)
+ sd30=1.0D0/(d_gmcrim*dgmm)
+ sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
+ sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
+ else
+ d30=0.0D0
+ d31=0.0D0
+ d32=0.0D0
+ d33=0.0D0
+ d34=0.0D0
+ d35=0.0D0
+ sd30=0.0D0
+ sd31=0.0D0
+ sd32=0.0D0
+ endif
+
+c write(*,*) 'inside angvectors4'
+
+cc**********************************************************************
+ atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
+ & -(ciy*gmz-ciz*gmy)*d31
+ & -(gmmy*rimmz-gmmz*rimmy)*d32
+ atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
+ & -(ciz*gmx-cix*gmz)*d31
+ & -(gmmz*rimmx-gmmx*rimmz)*d32
+ atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
+ & -(cix*gmy-ciy*gmx)*d31
+ & -(gmmx*rimmy-gmmy*rimmx)*d32
+cc**********************************************************************
+ astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
+ & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
+ & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
+ & -ciy*rimy*gmmx-rimz*gmx*rimmz)
+ & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
+ & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
+ & -sd32*(gmmy*rimmz-rimmy*gmmz)
+
+ astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
+ & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
+ & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
+ & +gmz*rimy*rimmz-rimz*ciz*gmmy)
+ & -sd31*(gmcrim_x*(cix*rimy-gmz)
+ & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
+ & -sd32*(-gmmx*rimmz+rimmx*gmmz)
+
+ astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
+ & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
+ & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
+ & +rimz*ciy*gmmy+rimz*gmx*rimmx)
+ & -sd31*(gmcrim_x*(cix*rimz+gmy)
+ & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
+ & -sd32*(gmmx*rimmy-rimmx*gmmy)
+c**********************************************************************
+ apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
+ & -(gmmy*rimmz-gmmz*rimmy)*d34
+ & +rix*d35
+ apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
+ & -(gmmz*rimmx-gmmx*rimmz)*d34
+ & +riy*d35
+ apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
+ & -(gmmx*rimmy-gmmy*rimmx)*d34
+ & +riz*d35
+ endif
+
+ if(i.ge.4) then
+ if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
+ & .and.drim3.gt.0.0D0
+ & .and.d_gmmcrimm3.gt.0.0D0) then
+ d40=1.0D0/(dg30*dgmm)
+ d41=g3gmm/(dg30*dgmm3)
+ d42=1.0D0/(dg30*drim)
+ d43=g3rim_/(dg30*drim3)
+ sd40=1.0D0/(dg30*d_gmmcrimm)
+ sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
+ else
+ d40=0.0D0
+ d41=0.0D0
+ d42=0.0D0
+ d43=0.0D0
+ sd40=0.0D0
+ sd41=0.0D0
+ endif
+ atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
+ & -(gmmy*rimmz-gmmz*rimmy)*d41
+ atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
+ & -(gmmz*rimmx-gmmx*rimmz)*d41
+ atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
+ & -(gmmx*rimmy-gmmy*rimmx)*d41
+cc**********************************************************************
+ astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
+ & -g3z*rimmz*rimmx+rimmy**2*g3x)
+ & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
+ & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
+
+ astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
+ & -rimmx*rimmy*g3x+rimmz**2*g3y)
+ & -sd41*(-gmmcrimm_x*rimmx*rimmy
+ & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmy)
+
+c & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
+
+ astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
+ & +g3z*rimmx**2-rimmz*rimmy*g3y)
+ & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
+ & +gmmcrimm_z*(rimmy**2+rimmx**2))
+c**********************************************************************
+ apm3x(i)=g3x*d42-rimx*d43
+ apm3y(i)=g3y*d42-rimy*d43
+ apm3z(i)=g3z*d42-rimz*d43
+ endif
+ enddo
+c*******************************************************************************
+
+c write(*,*) 'inside angvectors5'
+
+c do i=inb-2,inb
+ do i=1,0
+ rimx=bx(i)-bx(i-1)
+ rimy=by(i)-by(i-1)
+ rimz=bz(i)-bz(i-1)
+ rimmx=bx(i-1)-bx(i-2)
+ rimmy=by(i-1)-by(i-2)
+ rimmz=bz(i-1)-bz(i-2)
+ rim3x=bx(i-2)-bx(i-3)
+ rim3y=by(i-2)-by(i-3)
+ rim3z=bz(i-2)-bz(i-3)
+ gmmx=rimmy*rimz-rimmz*rimy
+ gmmy=rimmz*rimx-rimmx*rimz
+ gmmz=rimmx*rimy-rimmy*rimx
+ g3x=rim3y*rimmz-rim3z*rimmy
+ g3y=rim3z*rimmx-rim3x*rimmz
+ g3z=rim3x*rimmy-rim3y*rimmx
+
+ dg30=sqrt(g3x**2+g3y**2+g3z**2)
+ g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
+ dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
+ dgmm3=dgmm**3
+ drim=dis(i-1,i)
+ drimm=dis(i-2,i-1)
+ drim3=drim**3
+ g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
+cc**********************************************************************
+ gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
+ gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
+ gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
+ d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
+ d_gmmcrimm3=d_gmmcrimm**3
+ gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
+ & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
+
+ if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
+ & .and.drim3.gt.0.0D0
+ & .and.d_gmmcrimm3.gt.0.0D0) then
+ d40=1.0D0/(dg30*dgmm)
+ d41=g3gmm/(dg30*dgmm3)
+ d42=1.0D0/(dg30*drim)
+ d43=g3rim_/(dg30*drim3)
+ sd40=1.0D0/(dg30*d_gmmcrimm)
+ sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
+ else
+ d40=0.0D0
+ d41=0.0D0
+ d42=0.0D0
+ d43=0.0D0
+ sd40=0.0D0
+ sd41=0.0D0
+ endif
+ atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
+ & -(gmmy*rimmz-gmmz*rimmy)*d41
+ atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
+ & -(gmmz*rimmx-gmmx*rimmz)*d41
+ atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
+ & -(gmmx*rimmy-gmmy*rimmx)*d41
+cc**********************************************************************
+ astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
+ & -g3z*rimmz*rimmx+rimmy**2*g3x)
+ & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
+ & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
+
+ astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
+ & -rimmx*rimmy*g3x+rimmz**2*g3y)
+ & -sd41*(-gmmcrimm_x*rimmx*rimmy
+ & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
+
+ astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
+ & +g3z*rimmx**2-rimmz*rimmy*g3y)
+ & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
+ & +gmmcrimm_z*(rimmy**2+rimmx**2))
+cc**********************************************************************
+ apm3x(i)=g3x*d42-rimx*d43
+ apm3y(i)=g3y*d42-rimy*d43
+ apm3z(i)=g3z*d42-rimz*d43
+
+ if(i.le.inb-1) then
+ ip=i+1
+ rix=bx(ip)-bx(i)
+ riy=by(ip)-by(i)
+ riz=bz(ip)-bz(i)
+ cix=bx(ip)-bx(i-1)
+ ciy=by(ip)-by(i-1)
+ ciz=bz(ip)-bz(i-1)
+ gmx=rimy*riz-rimz*riy
+ gmy=rimz*rix-rimx*riz
+ gmz=rimx*riy-rimy*rix
+ dgm=sqrt(gmx**2+gmy**2+gmz**2)
+ dgm3=dgm**3
+ dri=dis(i,i+1)
+ dri3=dri**3
+ gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
+ gmmr=gmmx*rix+gmmy*riy+gmmz*riz
+ gmcrim_x=gmy*rimz-gmz*rimy
+ gmcrim_y=gmz*rimx-gmx*rimz
+ gmcrim_z=gmx*rimy-gmy*rimx
+ d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
+ d_gmcrim3=d_gmcrim**3
+ gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
+ & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
+
+ if(dgm3.gt.0.0D0.and.
+ & dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
+ & .and.d_gmcrim3.gt.0.0D0) then
+ d30=1.0D0/(dgm*dgmm)
+ d31=gmmgm/(dgm3*dgmm)
+ d32=gmmgm/(dgm*dgmm3)
+ d33=1.0D0/(dgmm*dri)
+ d34=gmmr/(dgmm3*dri)
+ d35=gmmr/(dgmm*dri3)
+ sd30=1.0D0/(d_gmcrim*dgmm)
+ sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
+ sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
+
+ else
+ d30=0.0D0
+ d31=0.0D0
+ d32=0.0D0
+ d33=0.0D0
+ d34=0.0D0
+ d35=0.0D0
+ sd30=0.0D0
+ sd31=0.0D0
+ sd32=0.0D0
+ endif
+cc**********************************************************************
+ atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
+ & -(ciy*gmz-ciz*gmy)*d31
+ & -(gmmy*rimmz-gmmz*rimmy)*d32
+ atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
+ & -(ciz*gmx-cix*gmz)*d31
+ & -(gmmz*rimmx-gmmx*rimmz)*d32
+ atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
+ & -(cix*gmy-ciy*gmx)*d31
+ & -(gmmx*rimmy-gmmy*rimmx)*d32
+cc**********************************************************************
+ astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
+ & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
+ & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
+ & -ciy*rimy*gmmx-rimz*gmx*rimmz)
+ & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
+ & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
+ & -sd32*(gmmy*rimmz-rimmy*gmmz)
+
+ astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
+ & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
+ & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
+ & +gmz*rimy*rimmz-rimz*ciz*gmmy)
+ & -sd31*(gmcrim_x*(cix*rimy-gmz)
+ & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
+ & -sd32*(-gmmx*rimmz+rimmx*gmmz)
+
+ astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
+ & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
+ & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
+ & +rimz*ciy*gmmy+rimz*gmx*rimmx)
+ & -sd31*(gmcrim_x*(cix*rimz+gmy)
+ & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
+ & -sd32*(gmmx*rimmy-rimmx*gmmy)
+cc**********************************************************************
+ apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
+ & -(gmmy*rimmz-gmmz*rimmy)*d34
+ & +rix*d35
+ apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
+ & -(gmmz*rimmx-gmmx*rimmz)*d34
+ & +riy*d35
+ apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
+ & -(gmmx*rimmy-gmmy*rimmx)*d34
+ & +riz*d35
+ endif
+
+c write(*,*) 'inside angvectors6'
+
+ if(i.eq.inb-2) then
+ ipp=i+2
+ ripx=bx(ipp)-bx(ip)
+ ripy=by(ipp)-by(ip)
+ ripz=bz(ipp)-bz(ip)
+ cipx=bx(ipp)-bx(i)
+ cipy=by(ipp)-by(i)
+ cipz=bz(ipp)-bz(i)
+ gx=riy*ripz-riz*ripy
+ gy=riz*ripx-rix*ripz
+ gz=rix*ripy-riy*ripx
+ ggm=gmx*gx+gmy*gy+gmz*gz
+ gmrp=gmx*ripx+gmy*ripy+gmz*ripz
+ dg=sqrt(gx**2+gy**2+gz**2)
+ dg3=dg**3
+ drip=dis(i+1,i+2)
+ gcr_x=gy*riz-gz*riy
+ gcr_y=gz*rix-gx*riz
+ gcr_z=gx*riy-gy*rix
+ d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
+ d_gcr3=d_gcr**3
+ gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
+ & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
+ if(dgm3.gt.0.0D0.and.
+ & dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
+ & ) then
+ d20=1.0D0/(dg*dgm)
+ d21=ggm/(dgm3*dg)
+ d22=ggm/(dgm*dg3)
+ d23=1.0D0/(dgm*drip)
+ d24=gmrp/(dgm3*drip)
+ sd20=1.0D0/(d_gcr*dgm)
+ sd21=gcr__gm/(d_gcr3*dgm)
+ sd22=gcr__gm/(d_gcr*dgm3)
+ else
+ d20=0.0D0
+ d21=0.0D0
+ d22=0.0D0
+ d23=0.0D0
+ d24=0.0D0
+ sd20=0.0D0
+ sd21=0.0D0
+ sd22=0.0D0
+ endif
+ atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
+ & -(ciy*gmz-ciz*gmy)*d21
+ & +(ripy*gz-ripz*gy)*d22
+ atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
+ & -(ciz*gmx-cix*gmz)*d21
+ & +(ripz*gx-ripx*gz)*d22
+ atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
+ & -(cix*gmy-ciy*gmx)*d21
+ & +(ripx*gy-ripy*gx)*d22
+cc**********************************************************************
+ astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
+ & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
+ & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
+ & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
+ & +gcr_z*(-ripz*rix+gy))
+ & -sd22*(-gmy*ciz+gmz*ciy)
+
+ astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
+ & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
+ & +riz*ripz*gmy)
+ & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
+ & -gcr_z*(ripz*riy+gx))
+ & -sd22*(gmx*ciz-gmz*cix)
+
+ astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
+ & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
+ & -riz*gx*cix)
+ & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
+ & +gcr_z*(ripy*riy+ripx*rix))
+ & -sd22*(-gmx*ciy+gmy*cix)
+cc**********************************************************************
+c
+ apmx(i)=(ciy*ripz-ripy*ciz)*d23
+ & -(ciy*gmz-ciz*gmy)*d24
+ apmy(i)=(ciz*ripx-ripz*cix)*d23
+ & -(ciz*gmx-cix*gmz)*d24
+ apmz(i)=(cix*ripy-ripx*ciy)*d23
+ & -(cix*gmy-ciy*gmx)*d24
+
+ endif
+ enddo
+
+ return
+ end
+c END of angvectors
+c-------------------------------------------------------------------------------
+C---------------------------------------------------------------------------------
+ subroutine sheetforce(nca,wshet)
+ implicit none
+C JYLEE
+c this should be matched with dfa.fcm
+ integer maxca
+ parameter(maxca=800)
+cc**********************************************************************
+ integer nca
+ integer i,k
+ integer inb,nmax,iselect
+
+c real*8 dfaexp(15001)
+
+ real*8 vbeta,vbetp,vbetm
+ real*8 shefx(maxca,12)
+ real*8 shefy(maxca,12),shefz(maxca,12)
+ real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
+ real*8 vbet(maxca,maxca)
+ real*8 wshet(maxca,maxca)
+ real*8 bx(maxca),by(maxca),bz(maxca)
+
+ common /sheca/ bx,by,bz
+ common /phys1/ inb,nmax,iselect
+ common /shef/ shefx,shefy,shefz
+ common /shee/ vbeta,vbet,vbetp,vbetm
+ common /shetf/ shetfx,shetfy,shetfz
+
+ inb=nca
+ do i=1,inb
+ shetfx(i)=0.0D0
+ shetfy(i)=0.0D0
+ shetfz(i)=0.0D0
+ enddo
+
+ do k=1,12
+ do i=1,inb
+ shefx(i,k)=0.0D0
+ shefy(i,k)=0.0D0
+ shefz(i,k)=0.0D0
+ enddo
+ enddo
+
+ call sheetene(nca,wshet)
+ call sheetforce1
+
+ 887 format(a,1x,i6,3x,f12.8)
+ 888 format(a,1x,i4,1x,i4,3x,f12.8)
+ 889 format(a,1x,i4,3x,f12.8)
+ !write(2,*) 'coord : '
+ do i=1,inb
+ !write(2,887) 'bx:',i,bx(i)
+ !write(2,887) 'by:',i,by(i)
+ !write(2,887) 'bz:',i,bz(i)
+ enddo
+ !write(2,*) 'After sheetforce1'
+ do i=1,inb
+ do k=1,12
+ !write(2,888) 'shefx :',i,k,shefx(i,k)
+ !write(2,888) 'shefy :',i,k,shefy(i,k)
+ !write(2,888) 'shefz :',i,k,shefz(i,k)
+ enddo
+ enddo
+
+ call sheetforce5
+
+ !write(2,*) 'After sheetforce5'
+ do i=1,inb
+ do k=1,12
+ !write(2,888) 'shefx :',i,k,shefx(i,k)
+ !write(2,888) 'shefy :',i,k,shefy(i,k)
+ !write(2,888) 'shefz :',i,k,shefz(i,k)
+ enddo
+ enddo
+
+ call sheetforce6
+
+ !write(2,*) 'After sheetforce6'
+ do i=1,inb
+ do k=1,12
+ !write(2,888) 'shefx :',i,k,shefx(i,k)
+ !write(2,888) 'shefy :',i,k,shefy(i,k)
+ !write(2,888) 'shefz :',i,k,shefz(i,k)
+ enddo
+ enddo
+
+ call sheetforce11
+
+ !write(2,*) 'After sheetforce11'
+ do i=1,inb
+ do k=1,12
+ !write(2,888) 'shefx :',i,k,shefx(i,k)
+ !write(2,888) 'shefy :',i,k,shefy(i,k)
+ !write(2,888) 'shefz :',i,k,shefz(i,k)
+ enddo
+ enddo
+
+ call sheetforce12
+
+ !write(2,*) 'After sheetforce12'
+ do i=1,inb
+ do k=1,12
+ !write(2,888) 'shefx :',i,k,shefx(i,k)
+ !write(2,888) 'shefy :',i,k,shefy(i,k)
+ !write(2,888) 'shefz :',i,k,shefz(i,k)
+ enddo
+ enddo
+
+ do i=1,inb
+ do k=1,12
+ shetfx(i)=shetfx(i)+shefx(i,k)
+ shetfy(i)=shetfy(i)+shefy(i,k)
+ shetfz(i)=shetfz(i)+shefz(i,k)
+ enddo
+ enddo
+ !write(2,*) 'Beta Finished'
+ do i=1,inb
+ !write(2,889) 'shetfx : ',i,shetfx(i)
+ !write(2,889) 'shetfy : ',i,shetfy(i)
+ !write(2,889) 'shetfz : ',i,shetfz(i)
+ enddo
+
+ return
+ end
+C end sheetforce
+c-------------------------------------------------------------------------------
+ subroutine sheetene(nca,wshet)
+ implicit none
+ integer maxca
+ parameter(maxca=800)
+cc******************************************************************************
+
+c real*8 dfaexp(15001)
+ real*8 dtmp1, dtmp2, dtmp3
+
+ real*8 vbet(maxca,maxca)
+ real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
+ real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
+ real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
+ real*8 pin1(maxca,maxca),pin2(maxca,maxca)
+ real*8 pin3(maxca,maxca),pin4(maxca,maxca)
+ real*8 pina1(maxca,maxca),pina2(maxca,maxca)
+ real*8 pina3(maxca,maxca),pina4(maxca,maxca)
+ real*8 cph(maxca),cth(maxca)
+ real*8 rx(maxca,maxca)
+ real*8 ry(maxca,maxca),rz(maxca,maxca)
+ real*8 bx(maxca),by(maxca),bz(maxca)
+ real*8 dis(maxca,maxca)
+ real*8 ulcos(maxca)
+cc**********************************************************************
+ real*8 astx(maxca),asty(maxca),astz(maxca)
+ real*8 astmx(maxca),astmy(maxca),astmz(maxca)
+ real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
+ real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
+ real*8 sth(maxca)
+ real*8 wshet(maxca,maxca)
+ real*8 dp45, dm45, w_beta
+ real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
+ integer nca
+ integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
+ real*8 uum, uup
+ real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
+
+ common /sheca/ bx,by,bz
+ common /phys1/ inb,nmax,iselect
+ common /kyori2/ dis
+ common /difvec/ rx,ry,rz
+ common /coscos/ cph,cth
+ common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ & c00,s00,ulnex,dnex
+ common /sheconst/ dp45,dm45,w_beta
+ common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
+ common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
+ common /shee/ vbeta,vbet,vbetp,vbetm
+ common /ulang/ ulcos
+cc**********************************************************************
+ common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
+ & astmmz,astm3x,astm3y,astm3z
+ common /sinsin/ sth
+
+ real*8 r_pair_mat(maxca,maxca)
+ci integer istrand(maxca,maxca)
+ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
+ci common /shetest/ istrand,istrand_p,istrand_m
+ common /beta_p/ r_pair_mat
+C-------------------------------------------------------------------------------
+ r_pair_mat = 0.0d0
+ do i=1,inb
+ do j=1,inb
+ r_pair_mat(i,j)=wshet(i,j)
+c write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
+ enddo
+ enddo
+c stop
+c
+ vbeta=0.0D0
+ vbetp=0.0D0
+ vbetm=0.0D0
+
+ do i=1,inb-7
+ do j=i+4,inb-3
+ ip=i+1
+ ipp=i+2
+ jp=j+1
+ jpp=j+2
+cc**********************************************************************
+ y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
+ & +(cth(j)*c00+sth(j)*s00-1.0D0)**2
+ y1=-0.5d0*y1/dca
+ y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
+ & +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
+ y2=-0.5d0*y2/dnex
+
+cdebug y2=0
+
+ y=y1+y2
+
+ci if(y.ge.-4) then
+ci istrand(i,j)=1
+ci else
+ci istrand(i,j)=0
+ci endif
+
+ci if(istrand(i,j).eq.1) then
+
+ yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
+ yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
+
+
+ pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
+ $ +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
+ pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
+ $ +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
+ pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
+ $ +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
+ pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
+ $ +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
+
+ yshe1=pin1(i,j)**2+pin2(i,j)**2
+ yshe1=-0.5d0*yshe1/dshe
+ yshe2=pin3(i,j)**2+pin4(i,j)**2
+ yshe2=-0.5d0*yshe2/dshe
+
+ci if((yshe1+yshe2).ge.-4) then
+ci istrand_p(i,j)=1
+ci else
+ci istrand_p(i,j)=0
+ci endif
+
+
+C write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
+C write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
+C write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
+C write(*,*) 'dis(i,j):',i,j,dis(i,j)
+C write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
+C write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
+C write(*,*) 'pin1:',pin1(i,j)
+C write(*,*) 'pin2:',pin2(i,j)
+C write(*,*) 'pin3:',pin3(i,j)
+C write(*,*) 'pin4:',pin4(i,j)
+
+C write(*,*) 'y:',y
+C write(*,*) 'yy1:',yy1
+C write(*,*) 'yy2:',yy2
+C write(*,*) 'yshe1:',yshe1
+C write(*,*) 'yshe2:',yshe2
+c
+
+ci if (istrand_p(i,j).eq.1) then
+
+cd yy1=0
+cd yy2=0
+cd yshe1=0
+cd yshe2=0
+ dtmp1 = y+yy1+yshe1
+ dtmp2 = y+yy2+yshe2
+ dtmp3 = y+yy1+yy2+yshe1+yshe2
+
+C write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
+C write(*,*)'2', y,yy1,yy2
+C write(*,*)'3', yshe1,yshe2
+
+cc if (dtmp3.le.-35.0d0) then
+c vbetap(i,j)=-dp45*exp(dtmp3)
+cc vbetap(i,j)=0.0d0
+cc else
+c vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
+ vbetap(i,j)=-dp45*exp(dtmp3)
+cc end if
+
+cc if (dtmp1.le.-35.0d0) then
+c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
+cc vbetap1(i,j)=0.0d0
+cc else
+c vbetap1(i,j)=-r_pair_mat(i+1,j+1)
+c $ *dfaexp(idint(-dtmp1*1000)+1)
+ vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
+cc end if
+
+cc if (dtmp2.le.-35.0d0) then
+C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
+cc vbetap2(i,j)=0.0d0
+cc else
+c vbetap2(i,j)=-r_pair_mat(i+2,j+2)
+c $ *dfaexp(idint(-dtmp2*1000)+1)
+ vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
+cc end if
+
+c vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
+c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
+c vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
+
+! write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
+! write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
+
+ci elseif (istrand_p(i,j).eq.0)then
+ci vbetap(i,j)=0
+ci vbetap1(i,j)=0
+ci vbetap2(i,j)=0
+ci endif
+
+ yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
+ yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
+
+ pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
+ $ +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
+ pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
+ $ +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
+ pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
+ $ +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
+ pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
+ $ +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
+
+ yshe1=pina1(i,j)**2+pina2(i,j)**2
+ yshe1=-0.5d0*yshe1/dshe
+ yshe2=pina3(i,j)**2+pina4(i,j)**2
+ yshe2=-0.5d0*yshe2/dshe
+
+ci if((yshe1+yshe2).ge.-4) then
+ci istrand_m(i,j)=1
+ci else
+ci istrand_m(i,j)=0
+ci endif
+
+
+C write(*,*) 'pina1:',pina1(i,j)
+C write(*,*) 'pina2:',pina2(i,j)
+C write(*,*) 'pina3:',pina3(i,j)
+C write(*,*) 'pina4:',pina4(i,j)
+C write(*,*) 'yshe1:',yshe1
+C write(*,*) 'yshe2:',yshe2
+C write(*,*) 'dshe:',dshe
+
+ci if (istrand_m(i,j).eq.1) then
+
+cd yy1=0
+cd yy2=0
+cd yshe1=0
+cd yshe2=0
+
+ dtmp3=y+yy1+yy2+yshe1+yshe2
+ dtmp1=y+yy1+yshe1
+ dtmp2=y+yy2+yshe2
+
+cc if(dtmp3 .le. -35.0d0) then
+c vbetam(i,j)=-dm45*exp(dtmp3)
+cc vbetam(i,j)=0.0d0
+cc else
+c vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
+ vbetam(i,j)=-dm45*exp(dtmp3)
+cc end if
+
+cc if(dtmp1 .le. -35.0d0) then
+c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
+cc vbetam1(i,j)=0.0d0
+cc else
+c vbetam1(i,j)=-r_pair_mat(i+1,j+2)
+c $ *dfaexp(idint(-dtmp1*1000)+1)
+ vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
+cc end if
+
+cc if(dtmp2.le.-35.0d0) then
+c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
+cc vbetam2(i,j)=0.0d0
+cc else
+c vbetam2(i,j)=-r_pair_mat(i+2,j+1)
+c $ *dfaexp(idint(-dtmp2*1000)+1)
+ vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
+cc end if
+
+ci elseif (istrand_m(i,j).eq.0)then
+ci vbetam(i,j)=0
+ci vbetam1(i,j)=0
+ci vbetam2(i,j)=0
+ci endif
+
+
+c vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
+c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
+c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
+
+! write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
+! write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
+
+ uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
+ uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
+
+c write(*,*) 'uup,uum:', uup, uum
+
+c uup=vbetap1(i,j)+vbetap2(i,j)
+c uum=vbetam1(i,j)+vbetam2(i,j)
+
+ vbet(i,j)=uup+uum
+ vbetp=vbetp+uup
+ vbetm=vbetm+uum
+ vbeta=vbeta+vbet(i,j)
+
+ci elseif(istrand(i,j).eq.0)then
+ci vbet(i,j)=0
+ci endif
+
+c write(*,*) 'uup,uum:',uup,uum
+c write(*,*) 'vbetap(i,j):',vbetap(i,j)
+c write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
+c write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
+c write(*,*) 'vbetam(i,j):',vbetam(i,j)
+c write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
+c write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
+c write(*,*) 'uup:',uup
+c write(*,*) 'uum:',uum
+c write(*,*) 'vbetp:',vbetp
+c write(*,*) 'vbetm:',vbetm
+c write(*,*) 'vbet(i,j):',vbet(i,j)
+c stop
+
+ enddo
+ enddo
+
+! do i=1,inb-7
+! do j=i+4,inb-3
+! write(*,*) 'I,J:', i,j
+! write(*,*) 'vbetap(i,j):',vbetap(i,j)
+! write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
+! write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
+! write(*,*) 'vbetam(i,j):',vbetam(i,j)
+! write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
+! write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
+! write(*,*) 'vbet(i,j):',vbet(i,j)
+! enddo
+! enddo
+
+ return
+ end
+c-------------------------------------------------------------------------------
+ subroutine sheetforce1
+ implicit none
+ integer maxca
+ parameter(maxca=800)
+cc**********************************************************************
+ real*8 vbet(maxca,maxca)
+ real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
+ real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
+ real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
+ real*8 cph(maxca),cth(maxca)
+ real*8 rx(maxca,maxca)
+ real*8 ry(maxca,maxca),rz(maxca,maxca)
+ real*8 bx(maxca),by(maxca),bz(maxca)
+ real*8 dis(maxca,maxca)
+ real*8 shefx(maxca,12)
+ real*8 shefy(maxca,12),shefz(maxca,12)
+ real*8 atx(maxca),aty(maxca),atz(maxca)
+ real*8 atmx(maxca),atmy(maxca),atmz(maxca)
+ real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
+ real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
+ real*8 apx(maxca),apy(maxca),apz(maxca)
+ real*8 apmx(maxca),apmy(maxca),apmz(maxca)
+ real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
+ real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
+ real*8 ulcos(maxca)
+ real*8 astx(maxca),asty(maxca),astz(maxca)
+ real*8 astmx(maxca),astmy(maxca),astmz(maxca)
+ real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
+ real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
+ real*8 sth(maxca)
+ real*8 w_beta,dp45, dm45
+ real*8 vbeta, vbetp, vbetm
+ real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ integer inb,nmax,iselect
+
+ common /phys1/ inb,nmax,iselect
+ common /kyori2/ dis
+ common /difvec/ rx,ry,rz
+ common /coscos/ cph,cth
+ common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ common /sheconst/ dp45,dm45,w_beta
+ common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
+ common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
+ $ atmmz,atm3x,atm3y,atm3z
+ common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
+ $ apmmz,apm3x,apm3y,apm3z
+ common /shef/ shefx,shefy,shefz
+ common /shee/ vbeta,vbet,vbetp,vbetm
+ common /ulang/ ulcos
+c c**********************************************************************
+ common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
+ $ astmmz,astm3x,astm3y,astm3z
+ common /sinsin/ sth
+C--------------------------------------------------------------------------------
+c local variables
+ integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
+ real*8 c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
+ real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
+ real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
+ real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
+C--------------------------------------------------------------------------------
+ do i=4,inb-4
+ im3=i-3
+ imm=i-2
+ im=i-1
+ c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
+ v1=0.0D0
+ do j=i+1,inb-3
+ v1=v1+vbet(im3,j)
+ enddo
+ cc1=(ulcos(imm)-ulnex)/dnex
+ dmm=cc1/(dis(imm,im)*dis(im,i))
+ dmm__=cc1*ulcos(imm)/dis(im,i)**2
+ fx=rx(imm,im)*dmm-rx(im,i)*dmm__
+ fy=ry(imm,im)*dmm-ry(im,i)*dmm__
+ fz=rz(imm,im)*dmm-rz(im,i)*dmm__
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
+ fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
+ fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
+ shefx(i,1)=fx*v1
+ shefy(i,1)=fy*v1
+ shefz(i,1)=fz*v1
+ enddo
+
+ do i=3,inb-5
+ imm=i-2
+ im=i-1
+ ip=i+1
+ c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
+ v2=0.0D0
+ do j=i+2,inb-3
+ v2=v2+vbet(imm,j)
+ enddo
+ cc1=(ulcos(imm)-ulnex)/dnex
+ cc2=(ulcos(im)-ulnex)/dnex
+ dmm1=cc1/(dis(imm,im)*dis(im,i))
+ dmm2=cc2/(dis(im,i)*dis(i,ip))
+ dmm1__=cc1*ulcos(imm)/dis(im,i)**2
+ dmm2_1=cc2*ulcos(im)/dis(im,i)**2
+ dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
+cc**********************************************************************
+ fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
+ $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
+ fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
+ $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
+ fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
+ $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
+ fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
+ fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
+ shefx(i,2)=fx*v2
+ shefy(i,2)=fy*v2
+ shefz(i,2)=fz*v2
+ enddo
+ do i=2,inb-6
+ im=i-1
+ ip=i+1
+ ipp=i+2
+ c3=(cth(im)*c00+sth(im)*s00-1)/dca
+ v3=0.0D0
+ do j=i+3,inb-3
+ v3=v3+vbet(im,j)
+ enddo
+ cc2=(ulcos(im)-ulnex)/dnex
+ cc3=(ulcos(i)-ulnex)/dnex
+ dmm2=cc2/(dis(im,i)*dis(i,ip))
+ dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
+ dmm2_1=cc2*ulcos(im)/dis(im,i)**2
+ dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
+ dmm3__=cc3*ulcos(i)/dis(i,ip)**2
+ fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
+ $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
+ fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
+ $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
+ fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
+ $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
+ fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
+ fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
+ shefx(i,3)=fx*v3
+ shefy(i,3)=fy*v3
+ shefz(i,3)=fz*v3
+ enddo
+ do i=1,inb-7
+ ip=i+1
+ ipp=i+2
+ c4=(cth(i)*c00+sth(i)*s00-1)/dca
+ v4=0.0D0
+ do j=i+4,inb-3
+ v4=v4+vbet(i,j)
+ enddo
+ cc3=(ulcos(i)-ulnex)/dnex
+ dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
+ dmm3__=cc3*ulcos(i)/dis(i,ip)**2
+ fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
+ fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
+ fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atx(i)*c00+astx(i)*s00)*c4
+ fy=fy+(aty(i)*c00+asty(i)*s00)*c4
+ fz=fz+(atz(i)*c00+astz(i)*s00)*c4
+ shefx(i,4)=fx*v4
+ shefy(i,4)=fy*v4
+ shefz(i,4)=fz*v4
+ enddo
+ do j=8,inb
+ jm3=j-3
+ jmm=j-2
+ jm=j-1
+ c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
+ v7=0.0D0
+ do i=1,j-7
+ v7=v7+vbet(i,jm3)
+ enddo
+ cc7=(ulcos(jmm)-ulnex)/dnex
+ dmm=cc7/(dis(jmm,jm)*dis(jm,j))
+ dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
+ fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
+ fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
+ fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
+ fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
+ fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
+ shefx(j,7)=fx*v7
+ shefy(j,7)=fy*v7
+ shefz(j,7)=fz*v7
+ enddo
+ do j=7,inb-1
+ jm=j-1
+ jmm=j-2
+ jp=j+1
+ c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
+ v8=0.0D0
+ do i=1,j-6
+ v8=v8+vbet(i,jmm)
+ enddo
+ cc7=(ulcos(jmm)-ulnex)/dnex
+ cc8=(ulcos(jm)-ulnex)/dnex
+ dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
+ dmm8=cc8/(dis(jm,j)*dis(j,jp))
+ dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
+ dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
+ dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
+ fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
+ $ -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
+ fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
+ $ -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
+ fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
+ $ -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
+ fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
+ fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
+ shefx(j,8)=fx*v8
+ shefy(j,8)=fy*v8
+ shefz(j,8)=fz*v8
+ enddo
+
+ do j=6,inb-2
+ jm=j-1
+ jp=j+1
+ jpp=j+2
+ c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
+ v9=0.0D0
+ do i=1,j-5
+ v9=v9+vbet(i,jm)
+ enddo
+ cc8=(ulcos(jm)-ulnex)/dnex
+ cc9=(ulcos(j)-ulnex)/dnex
+ dmm8=cc8/(dis(jm,j)*dis(j,jp))
+ dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
+ dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
+ dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
+ dmm9__=cc9*ulcos(j)/dis(j,jp)**2
+ fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
+ $ -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
+ fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
+ $ -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
+ fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
+ $ -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
+ fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
+ fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
+ shefx(j,9)=fx*v9
+ shefy(j,9)=fy*v9
+ shefz(j,9)=fz*v9
+ enddo
+
+ do j=5,inb-3
+ jp=j+1
+ jpp=j+2
+ c10=(cth(j)*c00+sth(j)*s00-1)/dca
+ v10=0.0D0
+ do i=1,j-4
+ v10=v10+vbet(i,j)
+ enddo
+ cc9=(ulcos(j)-ulnex)/dnex
+ dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
+ dmm9__=cc9*ulcos(j)/dis(j,jp)**2
+ fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
+ fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
+ fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
+cd fx=0
+cd fy=0
+cd fz=0
+ fx=fx+(atx(j)*c00+astx(j)*s00)*c10
+ fy=fy+(aty(j)*c00+asty(j)*s00)*c10
+ fz=fz+(atz(j)*c00+astz(j)*s00)*c10
+ shefx(j,10)=fx*v10
+ shefy(j,10)=fy*v10
+ shefz(j,10)=fz*v10
+ enddo
+
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine sheetforce5
+ implicit none
+ integer maxca
+ parameter(maxca=800)
+cc**********************************************************************
+ real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
+ real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
+ real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
+ real*8 pin1(maxca,maxca),pin2(maxca,maxca)
+ real*8 pin3(maxca,maxca),pin4(maxca,maxca)
+ real*8 pina1(maxca,maxca),pina2(maxca,maxca)
+ real*8 pina3(maxca,maxca),pina4(maxca,maxca)
+ real*8 rx(maxca,maxca)
+ real*8 ry(maxca,maxca),rz(maxca,maxca)
+ real*8 bx(maxca),by(maxca),bz(maxca)
+ real*8 dis(maxca,maxca)
+ real*8 shefx(maxca,12),shefy(maxca,12)
+ real*8 shefz(maxca,12)
+ real*8 dp45,dm45,w_beta
+ real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ integer inb,nmax,iselect
+cc**********************************************************************
+ common /phys1/ inb,nmax,iselect
+ common /kyori2/ dis
+ common /difvec/ rx,ry,rz
+ common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ common /sheconst/ dp45,dm45,w_beta
+ common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
+ common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
+ common /shef/ shefx,shefy,shefz
+ci integer istrand(maxca,maxca)
+ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
+ci common /shetest/ istrand,istrand_p,istrand_m
+c********************************************************************************
+c local variables
+ integer i,imm,im,jp,jpp,j
+ real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
+ real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
+ real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
+ real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
+ real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
+ real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
+c********************************************************************************
+ do i=3,inb-5
+ imm=i-2
+ im=i-1
+ do j=i+2,inb-3
+ jp=j+1
+ jpp=j+2
+
+ci if(istrand(imm,j).eq.1
+ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
+
+
+ yy1=-(dis(i,jpp)-ulhb)/dlhb
+ y1x=rx(jpp,i)/dis(i,jpp)
+ y1y=ry(jpp,i)/dis(i,jpp)
+ y1z=rz(jpp,i)/dis(i,jpp)
+ y11x=yy1*y1x
+ y11y=yy1*y1y
+ y11z=yy1*y1z
+
+ yy33=1.0D0/(dis(im,jp)*dis(im,i))
+ yyy3=pin1(imm,j)/(dis(im,i)**2)
+ yy3=-pin1(imm,j)/dshe
+ y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
+ y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
+ y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
+
+ yy44=1.0D0/(dis(i,jpp)*dis(im,i))
+ yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
+ yyy4b=pin3(imm,j)/(dis(im,i)**2)
+ yy4=-pin3(imm,j)/dshe
+ y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
+ $ -yyy4b*rx(im,i))*yy4
+ y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
+ $ -yyy4b*ry(im,i))*yy4
+ y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
+ $ -yyy4b*rz(im,i))*yy4
+
+
+ yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
+ yyy5=pin4(imm,j)/(dis(i,jpp)**2)
+ yy5=-pin4(imm,j)/dshe
+ y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
+ y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
+ y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
+
+ sx=y11x+y3x+y4x+y5x
+ sy=y11y+y3y+y4y+y5y
+ sz=y11z+y3z+y4z+y5z
+
+ sx1=y3x
+ sy1=y3y
+ sz1=y3z
+ sx2=y11x+y4x+y5x
+ sy2=y11y+y4y+y5y
+ sz2=y11z+y4z+y5z
+
+ shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
+ $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
+ shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
+ $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
+ shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
+ $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
+
+! shefx(i,5)=shefx(i,5)
+! $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
+! shefy(i,5)=shefy(i,5)
+! $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
+! shefz(i,5)=shefz(i,5)
+! $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
+
+ yy6=-(dis(i,jp)-uldhb)/dldhb
+ y6x=rx(jp,i)/dis(i,jp)
+ y6y=ry(jp,i)/dis(i,jp)
+ y6z=rz(jp,i)/dis(i,jp)
+ y66x=yy6*y6x
+ y66y=yy6*y6y
+ y66z=yy6*y6z
+
+ yy88=1.0D0/(dis(im,jpp)*dis(im,i))
+ yyy8=pina1(imm,j)/(dis(im,i)**2)
+ yy8=-pina1(imm,j)/dshe
+ y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
+ y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
+ y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
+
+ yy99=1.0D0/(dis(jp,i)*dis(im,i))
+ yyy9a=pina3(imm,j)/(dis(jp,i)**2)
+ yyy9b=pina3(imm,j)/(dis(im,i)**2)
+ yy9=-pina3(imm,j)/dshe
+ y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
+ $ -yyy9b*rx(im,i))*yy9
+ y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
+ $ -yyy9b*ry(im,i))*yy9
+ y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
+ $ -yyy9b*rz(im,i))*yy9
+
+ yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
+ yyy10=pina4(imm,j)/(dis(jp,i)**2)
+ yy10=-pina4(imm,j)/dshe
+ y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
+ y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
+ y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
+
+ sx=y66x+y8x+y9x+y10x
+ sy=y66y+y8y+y9y+y10y
+ sz=y66z+y8z+y9z+y10z
+
+ sx1=y8x
+ sy1=y8y
+ sz1=y8z
+ sx2=y66x+y9x+y10x
+ sy2=y66y+y9y+y10y
+ sz2=y66z+y9z+y10z
+
+ shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
+ $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
+ shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
+ $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
+ shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
+ $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
+
+! shefx(i,5)=shefx(i,5)
+! $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
+! shefy(i,5)=shefy(i,5)
+! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
+! shefz(i,5)=shefz(i,5)
+! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
+
+ci endif
+
+ enddo
+ enddo
+
+ return
+ end
+c--------------------------------------------------------------------------c
+ subroutine sheetforce6
+ implicit none
+ integer maxca
+ parameter(maxca=800)
+cc**********************************************************************
+ real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
+ real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
+ real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
+ real*8 pin1(maxca,maxca),pin2(maxca,maxca)
+ real*8 pin3(maxca,maxca),pin4(maxca,maxca)
+ real*8 pina1(maxca,maxca),pina2(maxca,maxca)
+ real*8 pina3(maxca,maxca),pina4(maxca,maxca)
+ real*8 rx(maxca,maxca)
+ real*8 ry(maxca,maxca),rz(maxca,maxca)
+ real*8 bx(maxca),by(maxca),bz(maxca)
+ real*8 dis(maxca,maxca)
+ real*8 shefx(maxca,12),shefy(maxca,12)
+ real*8 shefz(maxca,12)
+ real*8 dp45,dm45,w_beta
+ real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ integer inb,nmax,iselect
+cc**********************************************************************
+ common /phys1/ inb,nmax,iselect
+ common /kyori2/ dis
+ common /difvec/ rx,ry,rz
+ common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ common /sheconst/ dp45,dm45,w_beta
+ common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
+ common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
+ common /shef/ shefx,shefy,shefz
+ci integer istrand(maxca,maxca)
+ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
+ci common /shetest/ istrand,istrand_p,istrand_m
+cc**********************************************************************
+C local variables
+ integer i,imm,im,jp,jpp,j,ip
+ real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
+ real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
+ real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
+ real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
+ real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
+ real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
+C********************************************************************************
+ do i=2,inb-6
+ ip=i+1
+ im=i-1
+ do j=i+3,inb-3
+ jp=j+1
+ jpp=j+2
+
+ci if(istrand(im,j).eq.1
+ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
+
+
+ yy1=-(dis(i,jp)-ulhb)/dlhb
+ y1x=rx(jp,i)/dis(i,jp)
+ y1y=ry(jp,i)/dis(i,jp)
+ y1z=rz(jp,i)/dis(i,jp)
+ y11x=yy1*y1x
+ y11y=yy1*y1y
+ y11z=yy1*y1z
+
+ yy33=1.0D0/(dis(i,jp)*dis(i,ip))
+ yyy3a=pin1(im,j)/(dis(i,jp)**2)
+ yyy3b=pin1(im,j)/(dis(i,ip)**2)
+ yy3=-pin1(im,j)/dshe
+ y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
+ $ +yyy3b*rx(i,ip))*yy3
+ y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
+ $ +yyy3b*ry(i,ip))*yy3
+ y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
+ $ +yyy3b*rz(i,ip))*yy3
+
+ yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
+ yyy4=pin2(im,j)/(dis(i,jp)**2)
+ yy4=-pin2(im,j)/dshe
+ y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
+ y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
+ y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
+
+ yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
+ yyy5=pin3(im,j)/(dis(i,ip)**2)
+ yy5=-pin3(im,j)/dshe
+ y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
+ y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
+ y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
+
+ sx=y11x+y3x+y4x+y5x
+ sy=y11y+y3y+y4y+y5y
+ sz=y11z+y3z+y4z+y5z
+
+ sx1=y11x+y3x+y4x
+ sy1=y11y+y3y+y4y
+ sz1=y11z+y3z+y4z
+ sx2=y5x
+ sy2=y5y
+ sz2=y5z
+
+ shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
+ $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
+ shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
+ $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
+ shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
+ $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
+! shefx(i,6)=shefx(i,6)
+! $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
+! shefy(i,6)=shefy(i,6)
+! $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
+! shefz(i,6)=shefz(i,6)
+! $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
+
+ yy6=-(dis(jpp,i)-uldhb)/dldhb
+ y6x=rx(jpp,i)/dis(jpp,i)
+ y6y=ry(jpp,i)/dis(jpp,i)
+ y6z=rz(jpp,i)/dis(jpp,i)
+ y66x=yy6*y6x
+ y66y=yy6*y6y
+ y66z=yy6*y6z
+
+ yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
+ yyy8a=pina1(im,j)/(dis(i,jpp)**2)
+ yyy8b=pina1(im,j)/(dis(i,ip)**2)
+ yy8=-pina1(im,j)/dshe
+ y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
+ $ +yyy8b*rx(i,ip))*yy8
+ y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
+ $ +yyy8b*ry(i,ip))*yy8
+ y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
+ $ +yyy8b*rz(i,ip))*yy8
+
+ yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
+ yyy9=pina2(im,j)/(dis(i,jpp)**2)
+ yy9=-pina2(im,j)/dshe
+ y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
+ y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
+ y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
+
+ yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
+ yyy10=pina3(im,j)/(dis(i,ip)**2)
+ yy10=-pina3(im,j)/dshe
+ y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
+ y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
+ y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
+
+ sx=y66x+y8x+y9x+y10x
+ sy=y66y+y8y+y9y+y10y
+ sz=y66z+y8z+y9z+y10z
+
+ sx1=y66x+y8x+y9x
+ sy1=y66y+y8y+y9y
+ sz1=y66z+y8z+y9z
+ sx2=y10x
+ sy2=y10y
+ sz2=y10z
+
+ shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
+ $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
+ shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
+ $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
+ shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
+ $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
+
+! shefx(i,6)=shefx(i,6)
+! $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
+! shefy(i,6)=shefy(i,6)
+! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
+! shefz(i,6)=shefz(i,6)
+! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
+
+ci endif
+
+ enddo
+ enddo
+
+ return
+ end
+c-----------------------------------------------------------------------
+ subroutine sheetforce11
+ implicit none
+ integer maxca
+ parameter(maxca=800)
+cc**********************************************************************
+ real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
+ real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
+ real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
+ real*8 pin1(maxca,maxca),pin2(maxca,maxca)
+ real*8 pin3(maxca,maxca),pin4(maxca,maxca)
+ real*8 pina1(maxca,maxca),pina2(maxca,maxca)
+ real*8 pina3(maxca,maxca),pina4(maxca,maxca)
+ real*8 rx(maxca,maxca)
+ real*8 ry(maxca,maxca),rz(maxca,maxca)
+ real*8 bx(maxca),by(maxca),bz(maxca)
+ real*8 dis(maxca,maxca)
+ real*8 shefx(maxca,12),shefy(maxca,12)
+ real*8 shefz(maxca,12)
+ real*8 dp45,dm45,w_beta
+ real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ integer inb,nmax,iselect
+cc**********************************************************************
+ common /phys1/ inb,nmax,iselect
+ common /kyori2/ dis
+ common /difvec/ rx,ry,rz
+ common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ common /sheconst/ dp45,dm45,w_beta
+ common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
+ common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
+ common /shef/ shefx,shefy,shefz
+ci integer istrand(maxca,maxca)
+ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
+ci common /shetest/ istrand,istrand_p,istrand_m
+C********************************************************************************
+C local variables
+ integer j,jm,jmm,ip,i,ipp
+ real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
+ real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
+ real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
+ real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
+ real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
+ real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
+C********************************************************************************
+
+ do j=7,inb-1
+ jm=j-1
+ jmm=j-2
+ do i=1,j-6
+ ip=i+1
+ ipp=i+2
+
+ci if(istrand(i,jmm).eq.1
+ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
+
+
+ yy1=-(dis(ipp,j)-ulhb)/dlhb
+ y1x=rx(ipp,j)/dis(ipp,j)
+ y1y=ry(ipp,j)/dis(ipp,j)
+ y1z=rz(ipp,j)/dis(ipp,j)
+ y11x=yy1*y1x
+ y11y=yy1*y1y
+ y11z=yy1*y1z
+
+ yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
+ yyy3=pin2(i,jmm)/(dis(jm,j)**2)
+ yy3=-pin2(i,jmm)/dshe
+ y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
+ y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
+ y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
+
+ yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
+ yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
+ yy4=-pin3(i,jmm)/dshe
+ y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
+ y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
+ y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
+
+ yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
+ yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
+ yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
+ yy5=-pin4(i,jmm)/dshe
+ y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
+ $ -yyy5b*rx(jm,j))*yy5
+ y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
+ $ -yyy5b*ry(jm,j))*yy5
+ y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
+ $ -yyy5b*rz(jm,j))*yy5
+
+ sx=y11x+y3x+y4x+y5x
+ sy=y11y+y3y+y4y+y5y
+ sz=y11z+y3z+y4z+y5z
+
+ sx1=y3x
+ sy1=y3y
+ sz1=y3z
+ sx2=y11x+y4x+y5x
+ sy2=y11y+y4y+y5y
+ sz2=y11z+y4z+y5z
+
+ shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
+ $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
+ shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
+ $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
+ shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
+ $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
+
+! shefx(j,11)=shefx(j,11)
+! $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
+! shefy(j,11)=shefy(j,11)
+! $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
+! shefz(j,11)=shefz(j,11)
+! $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
+
+ yy6=-(dis(ip,j)-uldhb)/dldhb
+ y6x=rx(ip,j)/dis(ip,j)
+ y6y=ry(ip,j)/dis(ip,j)
+ y6z=rz(ip,j)/dis(ip,j)
+ y66x=yy6*y6x
+ y66y=yy6*y6y
+ y66z=yy6*y6z
+
+ yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
+ yyy8=pina1(i,jmm)/(dis(ip,j)**2)
+ yy8=-pina1(i,jmm)/dshe
+ y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
+ y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
+ y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
+
+ yy99=1.0D0/(dis(ip,j)*dis(jm,j))
+ yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
+ yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
+ yy9=-pina2(i,jmm)/dshe
+ y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
+ $ -yyy9b*rx(jm,j))*yy9
+ y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
+ $ -yyy9b*ry(jm,j))*yy9
+ y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
+ $ -yyy9b*rz(jm,j))*yy9
+
+ yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
+ yyy10=pina4(i,jmm)/(dis(jm,j)**2)
+ yy10=-pina4(i,jmm)/dshe
+ y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
+ y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
+ y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
+
+ sx=y66x+y8x+y9x+y10x
+ sy=y66y+y8y+y9y+y10y
+ sz=y66z+y8z+y9z+y10z
+
+ sx1=y66x+y8x+y9x
+ sy1=y66y+y8y+y9y
+ sz1=y66z+y8z+y9z
+ sx2=y10x
+ sy2=y10y
+ sz2=y10z
+
+ shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
+ $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
+ shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
+ $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
+ shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
+ $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
+
+! shefx(j,11)=shefx(j,11)
+! $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
+! shefy(j,11)=shefy(j,11)
+! $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
+! shefz(j,11)=shefz(j,11)
+! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
+
+ci endif
+
+ enddo
+ enddo
+
+ return
+ end
+c-----------------------------------------------------------------------
+ subroutine sheetforce12
+ implicit none
+ integer maxca
+ parameter(maxca=800)
+cc**********************************************************************
+ real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
+ real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
+ real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
+ real*8 pin1(maxca,maxca),pin2(maxca,maxca)
+ real*8 pin3(maxca,maxca),pin4(maxca,maxca)
+ real*8 pina1(maxca,maxca),pina2(maxca,maxca)
+ real*8 pina3(maxca,maxca),pina4(maxca,maxca)
+ real*8 rx(maxca,maxca)
+ real*8 ry(maxca,maxca),rz(maxca,maxca)
+ real*8 bx(maxca),by(maxca),bz(maxca)
+ real*8 dis(maxca,maxca)
+ real*8 shefx(maxca,12),shefy(maxca,12)
+ real*8 shefz(maxca,12)
+ real*8 dp45,dm45,w_beta
+ real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ integer inb,nmax,iselect
+cc**********************************************************************
+ common /phys1/ inb,nmax,iselect
+ common /kyori2/ dis
+ common /difvec/ rx,ry,rz
+ common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
+ $ c00,s00,ulnex,dnex
+ common /sheconst/ dp45,dm45,w_beta
+ common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
+ common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
+ common /shef/ shefx,shefy,shefz
+ci integer istrand(maxca,maxca)
+ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
+ci common /shetest/ istrand,istrand_p,istrand_m
+cc**********************************************************************
+C local variables
+ integer j,jm,jmm,ip,i,ipp,jp
+ real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
+ real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
+ real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
+ real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
+ real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
+!c*************************************************************************c
+ do j=6,inb-2
+ jp=j+1
+ jm=j-1
+ do i=1,j-5
+ ip=i+1
+ ipp=i+2
+
+ci if(istrand(i,jm).eq.1
+ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
+
+
+ yy1=-(dis(ip,j)-ulhb)/dlhb
+ y1x=rx(ip,j)/dis(ip,j)
+ y1y=ry(ip,j)/dis(ip,j)
+ y1z=rz(ip,j)/dis(ip,j)
+ y11x=y1x*yy1
+ y11y=y1y*yy1
+ y11z=y1z*yy1
+
+ yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
+ yyy3=pin1(i,jm)/(dis(ip,j)**2)
+ yy3=-pin1(i,jm)/dshe
+ y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
+ y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
+ y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
+ yy44=1.0D0/(dis(ip,j)*dis(j,jp))
+
+ yyy4a=pin2(i,jm)/(dis(ip,j)**2)
+ yyy4b=pin2(i,jm)/(dis(j,jp)**2)
+ yy4=-pin2(i,jm)/dshe
+ y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
+ $ +yyy4b*rx(j,jp))*yy4
+ y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
+ $ +yyy4b*ry(j,jp))*yy4
+ y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
+ $ +yyy4b*rz(j,jp))*yy4
+
+ yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
+ yyy5=pin4(i,jm)/(dis(j,jp)**2)
+ yy5=-pin4(i,jm)/dshe
+ y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
+ y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
+ y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
+
+ sx=y11x+y3x+y4x+y5x
+ sy=y11y+y3y+y4y+y5y
+ sz=y11z+y3z+y4z+y5z
+
+ sx1=y11x+y3x+y4x
+ sy1=y11y+y3y+y4y
+ sz1=y11z+y3z+y4z
+ sx2=y5x
+ sy2=y5y
+ sz2=y5z
+
+ shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
+ $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
+ shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
+ $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
+ shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
+ $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
+
+! shefx(j,12)=shefx(j,12)
+! $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
+! shefy(j,12)=shefy(j,12)
+! $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
+! shefz(j,12)=shefz(j,12)
+! $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
+
+ yy6=-(dis(ipp,j)-uldhb)/dldhb
+ y6x=rx(ipp,j)/dis(ipp,j)
+ y6y=ry(ipp,j)/dis(ipp,j)
+ y6z=rz(ipp,j)/dis(ipp,j)
+ y66x=yy6*y6x
+ y66y=yy6*y6y
+ y66z=yy6*y6z
+
+ yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
+ yyy8=pina2(i,jm)/(dis(j,jp)**2)
+ yy8=-pina2(i,jm)/dshe
+ y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
+ y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
+ y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
+
+ yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
+ yyy9=pina3(i,jm)/(dis(j,ipp)**2)
+ yy9=-pina3(i,jm)/dshe
+ y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
+ y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
+ y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
+
+ yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
+ yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
+ yyy10b=pina4(i,jm)/(dis(j,jp)**2)
+ yy10=-pina4(i,jm)/dshe
+ y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
+ $ +yyy10b*rx(j,jp))*yy10
+ y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
+ $ +yyy10b*ry(j,jp))*yy10
+ y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
+ $ +yyy10b*rz(j,jp))*yy10
+
+ sx=y66x+y8x+y9x+y10x
+ sy=y66y+y8y+y9y+y10y
+ sz=y66z+y8z+y9z+y10z
+
+ sx1=y8x
+ sy1=y8y
+ sz1=y8z
+ sx2=y66x+y9x+y10x
+ sy2=y66y+y9y+y10y
+ sz2=y66z+y9z+y10z
+
+ shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
+ $ -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
+ shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
+ $ -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
+ shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
+ $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
+
+ci endif
+
+ ENDDO
+ ENDDO
+
+ RETURN
+ END
+C===============================================================================
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+
+C BARTEK for dfa test!
+ if (wdfa_dist.gt.0) call edfad(edfadis)
+c print*, 'edfad is finished!', edfadis
+ if (wdfa_tor.gt.0) call edfat(edfator)
+c print*, 'edfat is finished!', edfator
+ if (wdfa_nei.gt.0) call edfan(edfanei)
+c print*, 'edfan is finished!', edfanei
+ if (wdfa_beta.gt.0) call edfab(edfabet)
+c print*, 'edfab is finished!', edfabet
+C stop
+C BARTEK
+
c print *,"Processor",myrank," computed USCSC"
#ifdef TIMING
#ifdef MPI
energia(21)=esccor
energia(22)=evdw_p
energia(23)=evdw_m
+ energia(24)=edfadis
+ energia(25)=edfator
+ energia(26)=edfanei
+ energia(27)=edfabet
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
c print *," Processor",myrank," left SUM_ENERGY"
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ edfadis=energia(24)
+ edfator=energia(25)
+ edfanei=energia(26)
+ edfabet=energia(27)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor
+ & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+ & +wdfa_beta*edfabet
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor
+ & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+ & +wdfa_beta*edfabet
+
#endif
energia(0)=etot
c detecting NaNQ
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
+
enddo
enddo
#else
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
+
enddo
enddo
#endif
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
+
+
enddo
enddo
#endif
& +wturn3*gel_loc_turn3(i)
& +wturn6*gel_loc_turn6(i)
& +wel_loc*gel_loc_loc(i)
+ & +wsccor*gsccor_loc(i)
enddo
#ifdef DEBUG
write (iout,*) "gloc after adding corr"
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+C Bartek
+ edfadis = energia(24)
+ edfator = energia(25)
+ edfanei = energia(26)
+ edfabet = energia(27)
+
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
& edihcnstr,ebr*nss,
- & Uconst,etot
+ & Uconst,edfadis,edfator,edfanei,edfabet,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
+ & 'EDFAD= ',1pE16.6,' (DFA distance energy)'/
+ & 'EDFAT= ',1pE16.6,' (DFA torsion energy)'/
+ & 'EDFAN= ',1pE16.6,' (DFA NCa energy)'/
+ & 'EDFAB= ',1pE16.6,' (DFA Beta energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ebr*nss,Uconst,etot
+ & ebr*nss,
+ & Uconst,edfadis,edfator,edfanei,edfabet,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
+ & 'EDFAD= ',1pE16.6,' (DFA distance energy)'/
+ & 'EDFAT= ',1pE16.6,' (DFA torsion energy)'/
+ & 'EDFAN= ',1pE16.6,' (DFA NCa energy)'/
+ & 'EDFAB= ',1pE16.6,' (DFA Beta energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
& "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
& "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
& "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB ","EVDWPP ",
- & "ESTR ","EVDW2_14 ","UCONST ", " ","ESCCOR"," "," "/
+ & "ESTR ","EVDW2_14 ","UCONST ", " ","ESCCOR"," "," ",
+ & "DFA DIS","DFA TOR","DFA NEI","DFA BET"/
data wname /
& "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
& "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
& "WSTRAIN","WVDWPP","WBOND","SCAL14"," "," ","WSCCOR",
- & " "," "/
- data nprint_ene /20/
+ & " "," ","WDFAD","WDFAT","WDFAN","WDFAB"/
+ data nprint_ene /24/
data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
- & 21,0,0,0/
+ & 21,24,25,26,27,0,0,0/
end
c---------------------------------------------------------------------------
subroutine init_int_table
CCC
C
read (isccor,*,end=113,err=113) nsccortyp
+#ifdef SCCORPDB
read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
do i=-ntyp,-1
isccortyp(i)=-isccortyp(-i)
call reada(weightcard,'WTORD',wtor_d,1.0D0)
call reada(weightcard,'WANG',wang,1.0D0)
call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+C Bartek
+ call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
+ call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
+ call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
+ call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
+C
call reada(weightcard,'SCAL14',scal14,0.4D0)
call reada(weightcard,'SCALSCP',scalscp,1.0d0)
call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
weights(18)=scal14
weights(21)=wsccor
endif
+C Bartek
+ weights(24)=wdfa_dist
+ weights(25)=wdfa_tor
+ weights(26)=wdfa_nei
+ weights(27)=wdfa_beta
if(me.eq.king.or..not.out1file)
& write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
- & wturn4,wturn6
+ & wturn4,wturn6,
+ & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
+
10 format (/'Energy-term weights (unscaled):'//
& 'WSCC= ',f10.6,' (SC-SC)'/
& 'WSCP= ',f10.6,' (SC-p)'/
& 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
& 'WTURN3= ',f10.6,' (turns, 3rd order)'/
& 'WTURN4= ',f10.6,' (turns, 4th order)'/
- & 'WTURN6= ',f10.6,' (turns, 6th order)')
+ & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+ & 'WDFA_D= ',f10.6,' (DFA, distance)' /
+ & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
+ & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
+ & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
+
if(me.eq.king.or..not.out1file)then
if (wcorr4.gt.0.0d0) then
write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
if(me.eq.king.or..not.out1file)
& write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
- & wturn4,wturn6
+ & wturn4,wturn6,
+ & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
+
22 format (/'Energy-term weights (scaled):'//
& 'WSCC= ',f10.6,' (SC-SC)'/
& 'WSCP= ',f10.6,' (SC-p)'/
& 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
& 'WTURN3= ',f10.6,' (turns, 3rd order)'/
& 'WTURN4= ',f10.6,' (turns, 4th order)'/
- & 'WTURN6= ',f10.6,' (turns, 6th order)')
+ & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+ & 'WDFA_D= ',f10.6,' (DFA, distance)' /
+ & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
+ & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
+ & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
+
if(me.eq.king.or..not.out1file)
& write (iout,*) "Reference temperature for weights calculation:",
& temp0
cd print *,'NNT=',NNT,' NCT=',NCT
if (itype(1).eq.ntyp1) nnt=2
if (itype(nres).eq.ntyp1) nct=nct-1
+
+C Juyong:READ init_vars
+C Initialize variables!
+C Juyong:READ read_info
+C READ fragment information!!
+C both routines should be in dfa.F file!!
+
+ if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
+ & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
+ call init_dfa_vars
+ print*, 'init_dfa_vars finished!'
+ call read_dfa_info
+ print*, 'read_dfa_info finished!'
+ endif
+C
+C
+
+
if (pdbref) then
if(me.eq.king.or..not.out1file)
& write (iout,'(a,i3)') 'nsup=',nsup
projects / unres.git / commitdiff