2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 include 'COMMON.SHIELD'
85 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
86 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
88 character*320 controlcard
93 read (INP,'(a)') titel
94 call card_concat(controlcard)
95 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
96 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
97 call reada(controlcard,'SEED',seed,0.0D0)
98 call random_init(seed)
99 C Set up the time limit (caution! The time must be input in minutes!)
100 read_cart=index(controlcard,'READ_CART').gt.0
101 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
102 C this variable with_theta_constr is the variable which allow to read and execute the
103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
104 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
105 write (iout,*) "with_theta_constr ",with_theta_constr
106 call readi(controlcard,'SYM',symetr,1)
107 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
108 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
109 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
110 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
111 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
112 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
113 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
114 call reada(controlcard,'DRMS',drms,0.1D0)
115 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
117 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
118 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
119 write (iout,'(a,f10.1)')'DRMS = ',drms
120 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
121 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
123 call readi(controlcard,'NZ_START',nz_start,0)
124 call readi(controlcard,'NZ_END',nz_end,0)
125 call readi(controlcard,'IZ_SC',iz_sc,0)
127 safety = 60.0d0*safety
130 call reada(controlcard,"T_BATH",t_bath,300.0d0)
131 minim=(index(controlcard,'MINIMIZE').gt.0)
132 dccart=(index(controlcard,'CART').gt.0)
133 overlapsc=(index(controlcard,'OVERLAP').gt.0)
134 overlapsc=.not.overlapsc
135 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
136 searchsc=.not.searchsc
137 sideadd=(index(controlcard,'SIDEADD').gt.0)
138 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
139 outpdb=(index(controlcard,'PDBOUT').gt.0)
140 outmol2=(index(controlcard,'MOL2OUT').gt.0)
141 pdbref=(index(controlcard,'PDBREF').gt.0)
142 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
143 indpdb=index(controlcard,'PDBSTART')
144 extconf=(index(controlcard,'EXTCONF').gt.0)
145 AFMlog=(index(controlcard,'AFM'))
146 selfguide=(index(controlcard,'SELFGUIDE'))
147 print *,'AFMlog',AFMlog,selfguide,"KUPA"
148 call readi(controlcard,'TUBEMOD',tubelog,0)
149 write (iout,*) TUBElog,"TUBEMODE"
150 call readi(controlcard,'IPRINT',iprint,0)
151 C SHIELD keyword sets if the shielding effect of side-chains is used
152 C 0 denots no shielding is used all peptide are equally despite the
153 C solvent accesible area
154 C 1 the newly introduced function
155 C 2 reseved for further possible developement
156 call readi(controlcard,'SHIELD',shield_mode,0)
157 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
158 write(iout,*) "shield_mode",shield_mode
160 call readi(controlcard,'TORMODE',tor_mode,0)
161 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
162 write(iout,*) "torsional and valence angle mode",tor_mode
163 call readi(controlcard,'MAXGEN',maxgen,10000)
164 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
165 call readi(controlcard,"KDIAG",kdiag,0)
166 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
167 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
168 & write (iout,*) "RESCALE_MODE",rescale_mode
169 split_ene=index(controlcard,'SPLIT_ENE').gt.0
170 if (index(controlcard,'REGULAR').gt.0.0D0) then
171 call reada(controlcard,'WEIDIS',weidis,0.1D0)
175 if (index(controlcard,'CHECKGRAD').gt.0) then
177 if (index(controlcard,'CART').gt.0) then
179 elseif (index(controlcard,'CARINT').gt.0) then
184 call reada(controlcard,'DELTA',aincr,1.0d-7)
185 elseif (index(controlcard,'THREAD').gt.0) then
187 call readi(controlcard,'THREAD',nthread,0)
188 if (nthread.gt.0) then
189 call reada(controlcard,'WEIDIS',weidis,0.1D0)
192 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
193 stop 'Error termination in Read_Control.'
195 else if (index(controlcard,'MCMA').gt.0) then
197 else if (index(controlcard,'MCEE').gt.0) then
199 else if (index(controlcard,'MULTCONF').gt.0) then
201 else if (index(controlcard,'MAP').gt.0) then
203 call readi(controlcard,'MAP',nmap,0)
204 else if (index(controlcard,'CSA').gt.0) then
206 crc else if (index(controlcard,'ZSCORE').gt.0) then
208 crc ZSCORE is rm from UNRES, modecalc=9 is available
211 cfcm else if (index(controlcard,'MCMF').gt.0) then
213 else if (index(controlcard,'SOFTREG').gt.0) then
215 else if (index(controlcard,'CHECK_BOND').gt.0) then
217 else if (index(controlcard,'TEST').gt.0) then
219 else if (index(controlcard,'MD').gt.0) then
221 else if (index(controlcard,'RE ').gt.0) then
225 lmuca=index(controlcard,'MUCA').gt.0
226 call readi(controlcard,'MUCADYN',mucadyn,0)
227 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
228 if (lmuca .and. (me.eq.king .or. .not.out1file ))
230 write (iout,*) 'MUCADYN=',mucadyn
231 write (iout,*) 'MUCASMOOTH=',muca_smooth
234 iscode=index(controlcard,'ONE_LETTER')
235 indphi=index(controlcard,'PHI')
236 indback=index(controlcard,'BACK')
237 iranconf=index(controlcard,'RAND_CONF')
238 i2ndstr=index(controlcard,'USE_SEC_PRED')
239 gradout=index(controlcard,'GRADOUT').gt.0
240 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
241 C DISTCHAINMAX become obsolete for periodic boundry condition
242 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
243 C Reading the dimensions of box in x,y,z coordinates
244 call reada(controlcard,'BOXX',boxxsize,100.0d0)
245 call reada(controlcard,'BOXY',boxysize,100.0d0)
246 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
247 c Cutoff range for interactions
248 call reada(controlcard,"R_CUT",r_cut,15.0d0)
249 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
250 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
251 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
252 if (lipthick.gt.0.0d0) then
253 bordliptop=(boxzsize+lipthick)/2.0
254 bordlipbot=bordliptop-lipthick
256 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
257 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
258 buflipbot=bordlipbot+lipbufthick
259 bufliptop=bordliptop-lipbufthick
260 if ((lipbufthick*2.0d0).gt.lipthick)
261 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
263 write(iout,*) "bordliptop=",bordliptop
264 write(iout,*) "bordlipbot=",bordlipbot
265 write(iout,*) "bufliptop=",bufliptop
266 write(iout,*) "buflipbot=",buflipbot
267 write (iout,*) "SHIELD MODE",shield_mode
268 if (TUBElog.gt.0) then
269 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
270 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
271 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
272 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
273 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
274 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
275 buftubebot=bordtubebot+tubebufthick
276 buftubetop=bordtubetop-tubebufthick
278 if (shield_mode.gt.0) then
280 C VSolvSphere the volume of solving sphere
282 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
283 C there will be no distinction between proline peptide group and normal peptide
284 C group in case of shielding parameters
285 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
286 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
287 write (iout,*) VSolvSphere,VSolvSphere_div
288 C long axis of side chain
290 long_r_sidechain(i)=vbldsc0(1,i)
291 short_r_sidechain(i)=sigma0(i)
295 if (me.eq.king .or. .not.out1file )
296 & write (iout,*) "DISTCHAINMAX",distchainmax
298 if(me.eq.king.or..not.out1file)
299 & write (iout,'(2a)') diagmeth(kdiag),
300 & ' routine used to diagonalize matrices.'
303 c--------------------------------------------------------------------------
304 subroutine read_REMDpar
308 implicit real*8 (a-h,o-z)
310 include 'COMMON.IOUNITS'
311 include 'COMMON.TIME1'
314 include 'COMMON.LANGEVIN'
316 include 'COMMON.LANGEVIN.lang0'
318 include 'COMMON.INTERACT'
319 include 'COMMON.NAMES'
321 include 'COMMON.REMD'
322 include 'COMMON.CONTROL'
323 include 'COMMON.SETUP'
325 character*320 controlcard
326 character*3200 controlcard1
327 integer iremd_m_total
329 if(me.eq.king.or..not.out1file)
330 & write (iout,*) "REMD setup"
332 call card_concat(controlcard)
333 call readi(controlcard,"NREP",nrep,3)
334 call readi(controlcard,"NSTEX",nstex,1000)
335 call reada(controlcard,"RETMIN",retmin,10.0d0)
336 call reada(controlcard,"RETMAX",retmax,1000.0d0)
337 mremdsync=(index(controlcard,'SYNC').gt.0)
338 call readi(controlcard,"NSYN",i_sync_step,100)
339 restart1file=(index(controlcard,'REST1FILE').gt.0)
340 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
341 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
342 if(max_cache_traj_use.gt.max_cache_traj)
343 & max_cache_traj_use=max_cache_traj
344 if(me.eq.king.or..not.out1file) then
345 cd if (traj1file) then
346 crc caching is in testing - NTWX is not ignored
347 cd write (iout,*) "NTWX value is ignored"
348 cd write (iout,*) " trajectory is stored to one file by master"
349 cd write (iout,*) " before exchange at NSTEX intervals"
351 write (iout,*) "NREP= ",nrep
352 write (iout,*) "NSTEX= ",nstex
353 write (iout,*) "SYNC= ",mremdsync
354 write (iout,*) "NSYN= ",i_sync_step
355 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
358 if (index(controlcard,'TLIST').gt.0) then
360 call card_concat(controlcard1)
361 read(controlcard1,*) (remd_t(i),i=1,nrep)
362 if(me.eq.king.or..not.out1file)
363 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
366 if (index(controlcard,'MLIST').gt.0) then
368 call card_concat(controlcard1)
369 read(controlcard1,*) (remd_m(i),i=1,nrep)
370 if(me.eq.king.or..not.out1file) then
371 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
374 iremd_m_total=iremd_m_total+remd_m(i)
376 write (iout,*) 'Total number of replicas ',iremd_m_total
379 if(me.eq.king.or..not.out1file)
380 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
383 c--------------------------------------------------------------------------
384 subroutine read_MDpar
388 implicit real*8 (a-h,o-z)
390 include 'COMMON.IOUNITS'
391 include 'COMMON.TIME1'
394 include 'COMMON.LANGEVIN'
396 include 'COMMON.LANGEVIN.lang0'
398 include 'COMMON.INTERACT'
399 include 'COMMON.NAMES'
401 include 'COMMON.SETUP'
402 include 'COMMON.CONTROL'
403 include 'COMMON.SPLITELE'
405 character*320 controlcard
407 call card_concat(controlcard)
408 call readi(controlcard,"NSTEP",n_timestep,1000000)
409 call readi(controlcard,"NTWE",ntwe,100)
410 call readi(controlcard,"NTWX",ntwx,1000)
411 call reada(controlcard,"DT",d_time,1.0d-1)
412 call reada(controlcard,"DVMAX",dvmax,2.0d1)
413 call reada(controlcard,"DAMAX",damax,1.0d1)
414 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
415 call readi(controlcard,"LANG",lang,0)
416 RESPA = index(controlcard,"RESPA") .gt. 0
417 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
418 ntime_split0=ntime_split
419 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
420 ntime_split0=ntime_split
421 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
422 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
423 rest = index(controlcard,"REST").gt.0
424 tbf = index(controlcard,"TBF").gt.0
425 usampl = index(controlcard,"USAMPL").gt.0
427 mdpdb = index(controlcard,"MDPDB").gt.0
428 call reada(controlcard,"T_BATH",t_bath,300.0d0)
429 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
430 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
431 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
432 if (count_reset_moment.eq.0) count_reset_moment=1000000000
433 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
434 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
435 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
436 if (count_reset_vel.eq.0) count_reset_vel=1000000000
437 large = index(controlcard,"LARGE").gt.0
438 print_compon = index(controlcard,"PRINT_COMPON").gt.0
439 rattle = index(controlcard,"RATTLE").gt.0
440 c if performing umbrella sampling, fragments constrained are read from the fragment file
446 if(me.eq.king.or..not.out1file) then
448 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
450 write (iout,'(a)') "The units are:"
451 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
452 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
453 & " acceleration: angstrom/(48.9 fs)**2"
454 write (iout,'(a)') "energy: kcal/mol, temperature: K"
456 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
457 write (iout,'(a60,f10.5,a)')
458 & "Initial time step of numerical integration:",d_time,
460 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
462 write (iout,'(2a,i4,a)')
463 & "A-MTS algorithm used; initial time step for fast-varying",
464 & " short-range forces split into",ntime_split," steps."
465 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
466 & r_cut," lambda",rlamb
468 write (iout,'(2a,f10.5)')
469 & "Maximum acceleration threshold to reduce the time step",
470 & "/increase split number:",damax
471 write (iout,'(2a,f10.5)')
472 & "Maximum predicted energy drift to reduce the timestep",
473 & "/increase split number:",edriftmax
474 write (iout,'(a60,f10.5)')
475 & "Maximum velocity threshold to reduce velocities:",dvmax
476 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
477 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
478 if (rattle) write (iout,'(a60)')
479 & "Rattle algorithm used to constrain the virtual bonds"
483 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
484 call reada(controlcard,"RWAT",rwat,1.4d0)
485 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
486 surfarea=index(controlcard,"SURFAREA").gt.0
487 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
488 if(me.eq.king.or..not.out1file)then
489 write (iout,'(/a,$)') "Langevin dynamics calculation"
492 & " with direct integration of Langevin equations"
493 else if (lang.eq.2) then
494 write (iout,'(a/)') " with TINKER stochasic MD integrator"
495 else if (lang.eq.3) then
496 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
497 else if (lang.eq.4) then
498 write (iout,'(a/)') " in overdamped mode"
500 write (iout,'(//a,i5)')
501 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
504 write (iout,'(a60,f10.5)') "Temperature:",t_bath
505 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
506 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
507 write (iout,'(a60,f10.5)')
508 & "Scaling factor of the friction forces:",scal_fric
509 if (surfarea) write (iout,'(2a,i10,a)')
510 & "Friction coefficients will be scaled by solvent-accessible",
511 & " surface area every",reset_fricmat," steps."
513 c Calculate friction coefficients and bounds of stochastic forces
514 eta=6*pi*cPoise*etawat
515 if(me.eq.king.or..not.out1file)
516 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
518 gamp=scal_fric*(pstok+rwat)*eta
519 stdfp=dsqrt(2*Rb*t_bath/d_time)
521 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
522 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
524 if(me.eq.king.or..not.out1file)then
525 write (iout,'(/2a/)')
526 & "Radii of site types and friction coefficients and std's of",
527 & " stochastic forces of fully exposed sites"
528 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
530 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
531 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
535 if(me.eq.king.or..not.out1file)then
536 write (iout,'(a)') "Berendsen bath calculation"
537 write (iout,'(a60,f10.5)') "Temperature:",t_bath
538 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
540 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
541 & count_reset_moment," steps"
543 & write (iout,'(a,i10,a)')
544 & "Velocities will be reset at random every",count_reset_vel,
548 if(me.eq.king.or..not.out1file)
549 & write (iout,'(a31)') "Microcanonical mode calculation"
551 if(me.eq.king.or..not.out1file)then
552 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
554 write(iout,*) "MD running with constraints."
555 write(iout,*) "Equilibration time ", eq_time, " mtus."
556 write(iout,*) "Constraining ", nfrag," fragments."
557 write(iout,*) "Length of each fragment, weight and q0:"
559 write (iout,*) "Set of restraints #",iset
561 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
562 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
564 write(iout,*) "constraints between ", npair, "fragments."
565 write(iout,*) "constraint pairs, weights and q0:"
567 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
568 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
570 write(iout,*) "angle constraints within ", nfrag_back,
571 & "backbone fragments."
572 write(iout,*) "fragment, weights:"
574 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
575 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
576 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
579 iset=mod(kolor,nset)+1
582 if(me.eq.king.or..not.out1file)
583 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
586 c------------------------------------------------------------------------------
589 C Read molecular data.
591 implicit real*8 (a-h,o-z)
597 include 'COMMON.IOUNITS'
600 include 'COMMON.INTERACT'
601 include 'COMMON.LOCAL'
602 include 'COMMON.NAMES'
603 include 'COMMON.CHAIN'
604 include 'COMMON.FFIELD'
605 include 'COMMON.SBRIDGE'
606 include 'COMMON.HEADER'
607 include 'COMMON.CONTROL'
608 include 'COMMON.DBASE'
609 include 'COMMON.THREAD'
610 include 'COMMON.CONTACTS'
611 include 'COMMON.TORCNSTR'
612 include 'COMMON.TIME1'
613 include 'COMMON.BOUNDS'
615 include 'COMMON.SETUP'
616 include 'COMMON.SHIELD'
617 character*4 sequence(maxres)
619 double precision x(maxvar)
620 character*256 pdbfile
621 character*400 weightcard
622 character*80 weightcard_t,ucase
623 dimension itype_pdb(maxres)
624 common /pizda/ itype_pdb
625 logical seq_comp,fail
626 double precision energia(0:n_ene)
632 C Read weights of the subsequent energy terms.
633 call card_concat(weightcard)
634 call reada(weightcard,'WLONG',wlong,1.0D0)
635 call reada(weightcard,'WSC',wsc,wlong)
636 call reada(weightcard,'WSCP',wscp,wlong)
637 call reada(weightcard,'WELEC',welec,1.0D0)
638 call reada(weightcard,'WVDWPP',wvdwpp,welec)
639 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
640 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
641 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
642 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
643 call reada(weightcard,'WTURN3',wturn3,1.0D0)
644 call reada(weightcard,'WTURN4',wturn4,1.0D0)
645 call reada(weightcard,'WTURN6',wturn6,1.0D0)
646 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
647 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
648 call reada(weightcard,'WBOND',wbond,1.0D0)
649 call reada(weightcard,'WTOR',wtor,1.0D0)
650 call reada(weightcard,'WTORD',wtor_d,1.0D0)
651 call reada(weightcard,'WANG',wang,1.0D0)
652 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
653 call reada(weightcard,'SCAL14',scal14,0.4D0)
654 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
655 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
656 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
657 call reada(weightcard,'TEMP0',temp0,300.0d0)
658 call reada(weightcard,'WSHIELD',wshield,1.0d0)
659 call reada(weightcard,'WLT',wliptran,0.0D0)
660 call reada(weightcard,'WTUBE',wtube,1.0D0)
661 if (index(weightcard,'SOFT').gt.0) ipot=6
662 C 12/1/95 Added weight for the multi-body term WCORR
663 call reada(weightcard,'WCORRH',wcorr,1.0D0)
664 if (wcorr4.gt.0.0d0) wcorr=wcorr4
684 if(me.eq.king.or..not.out1file)
685 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
686 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
688 10 format (/'Energy-term weights (unscaled):'//
689 & 'WSCC= ',f10.6,' (SC-SC)'/
690 & 'WSCP= ',f10.6,' (SC-p)'/
691 & 'WELEC= ',f10.6,' (p-p electr)'/
692 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
693 & 'WBOND= ',f10.6,' (stretching)'/
694 & 'WANG= ',f10.6,' (bending)'/
695 & 'WSCLOC= ',f10.6,' (SC local)'/
696 & 'WTOR= ',f10.6,' (torsional)'/
697 & 'WTORD= ',f10.6,' (double torsional)'/
698 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
699 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
700 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
701 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
702 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
703 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
704 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
705 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
706 & 'WTURN6= ',f10.6,' (turns, 6th order)')
707 if(me.eq.king.or..not.out1file)then
708 if (wcorr4.gt.0.0d0) then
709 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
710 & 'between contact pairs of peptide groups'
711 write (iout,'(2(a,f5.3/))')
712 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
713 & 'Range of quenching the correlation terms:',2*delt_corr
714 else if (wcorr.gt.0.0d0) then
715 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
716 & 'between contact pairs of peptide groups'
718 write (iout,'(a,f8.3)')
719 & 'Scaling factor of 1,4 SC-p interactions:',scal14
720 write (iout,'(a,f8.3)')
721 & 'General scaling factor of SC-p interactions:',scalscp
723 r0_corr=cutoff_corr-delt_corr
725 aad(i,1)=scalscp*aad(i,1)
726 aad(i,2)=scalscp*aad(i,2)
727 bad(i,1)=scalscp*bad(i,1)
728 bad(i,2)=scalscp*bad(i,2)
730 call rescale_weights(t_bath)
731 if(me.eq.king.or..not.out1file)
732 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
733 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
735 22 format (/'Energy-term weights (scaled):'//
736 & 'WSCC= ',f10.6,' (SC-SC)'/
737 & 'WSCP= ',f10.6,' (SC-p)'/
738 & 'WELEC= ',f10.6,' (p-p electr)'/
739 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
740 & 'WBOND= ',f10.6,' (stretching)'/
741 & 'WANG= ',f10.6,' (bending)'/
742 & 'WSCLOC= ',f10.6,' (SC local)'/
743 & 'WTOR= ',f10.6,' (torsional)'/
744 & 'WTORD= ',f10.6,' (double torsional)'/
745 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
746 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
747 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
748 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
749 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
750 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
751 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
752 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
753 & 'WTURN6= ',f10.6,' (turns, 6th order)')
754 if(me.eq.king.or..not.out1file)
755 & write (iout,*) "Reference temperature for weights calculation:",
757 call reada(weightcard,"D0CM",d0cm,3.78d0)
758 call reada(weightcard,"AKCM",akcm,15.1d0)
759 call reada(weightcard,"AKTH",akth,11.0d0)
760 call reada(weightcard,"AKCT",akct,12.0d0)
761 call reada(weightcard,"V1SS",v1ss,-1.08d0)
762 call reada(weightcard,"V2SS",v2ss,7.61d0)
763 call reada(weightcard,"V3SS",v3ss,13.7d0)
764 call reada(weightcard,"EBR",ebr,-5.50D0)
765 call reada(weightcard,"ATRISS",atriss,0.301D0)
766 call reada(weightcard,"BTRISS",btriss,0.021D0)
767 call reada(weightcard,"CTRISS",ctriss,1.001D0)
768 call reada(weightcard,"DTRISS",dtriss,1.001D0)
769 call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
770 write (iout,*) "ATRISS=", atriss
771 write (iout,*) "BTRISS=", btriss
772 write (iout,*) "CTRISS=", ctriss
773 write (iout,*) "DTRISS=", dtriss
774 if (shield_mode.gt.0) then
776 write (iout,*) "liscale not used in shield mode"
778 write (iout,*) "lipid scaling factor", lipscale
779 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
781 dyn_ss_mask(i)=.false.
785 dyn_ssbond_ij(i,j)=1.0d300
788 call reada(weightcard,"HT",Ht,0.0D0)
790 ss_depth=ebr/wsc-0.25*eps(1,1)
791 Ht=Ht/wsc-0.25*eps(1,1)
792 akcm=akcm*wstrain/wsc
793 akth=akth*wstrain/wsc
794 akct=akct*wstrain/wsc
795 v1ss=v1ss*wstrain/wsc
796 v2ss=v2ss*wstrain/wsc
797 v3ss=v3ss*wstrain/wsc
799 if (wstrain.ne.0.0) then
800 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
805 write (iout,*) "wshield,", wshield
806 if(me.eq.king.or..not.out1file) then
807 write (iout,*) "Parameters of the SS-bond potential:"
808 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
810 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
811 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
812 write (iout,*)" HT",Ht
813 print *,'indpdb=',indpdb,' pdbref=',pdbref
815 if (indpdb.gt.0 .or. pdbref) then
816 read(inp,'(a)') pdbfile
817 if(me.eq.king.or..not.out1file)
818 & write (iout,'(2a)') 'PDB data will be read from file ',
819 & pdbfile(:ilen(pdbfile))
820 open(ipdbin,file=pdbfile,status='old',err=33)
822 33 write (iout,'(a)') 'Error opening PDB file.'
825 c write (iout,*) 'Begin reading pdb data'
828 c write (iout,*) 'Finished reading pdb data'
830 if(me.eq.king.or..not.out1file)
831 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
832 & ' nstart_sup=',nstart_sup
834 itype_pdb(i)=itype(i)
838 nct=nstart_sup+nsup-1
839 call contact(.false.,ncont_ref,icont_ref,co)
842 if(me.eq.king.or..not.out1file)
843 & write(iout,*)'Adding sidechains'
847 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
850 do while (fail.and.nsi.le.maxsi)
851 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
854 if(fail) write(iout,*)'Adding sidechain failed for res ',
855 & i,' after ',nsi,' trials'
860 if (indpdb.eq.0) then
861 C Read sequence if not taken from the pdb file.
863 c print *,'nres=',nres
864 if (iscode.gt.0) then
865 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
867 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
869 C Convert sequence to numeric code
871 itype(i)=rescode(i,sequence(i),iscode)
873 C Assign initial virtual bond lengths
879 vbld(i+nres)=dsc(iabs(itype(i)))
880 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
881 c write (iout,*) "i",i," itype",itype(i),
882 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
886 c print '(20i4)',(itype(i),i=1,nres)
889 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
891 if (itype(i).eq.ntyp1) then
895 else if (iabs(itype(i+1)).ne.20) then
897 else if (iabs(itype(i)).ne.20) then
904 if(me.eq.king.or..not.out1file)then
905 write (iout,*) "ITEL"
907 write (iout,*) i,itype(i),itel(i)
909 print *,'Call Read_Bridge.'
912 C 8/13/98 Set limits to generating the dihedral angles
917 read (inp,*) ndih_constr
918 if (ndih_constr.gt.0) then
920 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
922 if(me.eq.king.or..not.out1file)then
924 & 'There are',ndih_constr,' constraints on phi angles.'
926 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
931 phi0(i)=deg2rad*phi0(i)
932 drange(i)=deg2rad*drange(i)
934 C if(me.eq.king.or..not.out1file)
935 C & write (iout,*) 'FTORS',ftors
938 phibound(1,ii) = phi0(i)-drange(i)
939 phibound(2,ii) = phi0(i)+drange(i)
942 C first setting the theta boundaries to 0 to pi
943 C this mean that there is no energy penalty for any angle occuring this can be applied
944 C for generate random conformation but is not implemented in this way
949 C begin reading theta constrains this is quartic constrains allowing to
950 C have smooth second derivative
951 if (with_theta_constr) then
952 C with_theta_constr is keyword allowing for occurance of theta constrains
953 read (inp,*) ntheta_constr
954 C ntheta_constr is the number of theta constrains
955 if (ntheta_constr.gt.0) then
957 read (inp,*) (itheta_constr(i),theta_constr0(i),
958 & theta_drange(i),for_thet_constr(i),
960 C the above code reads from 1 to ntheta_constr
961 C itheta_constr(i) residue i for which is theta_constr
962 C theta_constr0 the global minimum value
963 C theta_drange is range for which there is no energy penalty
964 C for_thet_constr is the force constant for quartic energy penalty
966 if(me.eq.king.or..not.out1file)then
968 & 'There are',ntheta_constr,' constraints on phi angles.'
970 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
976 theta_constr0(i)=deg2rad*theta_constr0(i)
977 theta_drange(i)=deg2rad*theta_drange(i)
979 C if(me.eq.king.or..not.out1file)
980 C & write (iout,*) 'FTORS',ftors
981 C do i=1,ntheta_constr
982 C ii = itheta_constr(i)
983 C thetabound(1,ii) = phi0(i)-drange(i)
984 C thetabound(2,ii) = phi0(i)+drange(i)
986 endif ! ntheta_constr.gt.0
987 endif! with_theta_constr
989 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
990 C write (iout,*) "with_dihed_constr ",with_dihed_constr
995 write (iout,'(a)') 'Boundaries in phi angle sampling:'
997 write (iout,'(a3,i5,2f10.1)')
998 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1004 cd print *,'NNT=',NNT,' NCT=',NCT
1005 if (itype(1).eq.ntyp1) nnt=2
1006 if (itype(nres).eq.ntyp1) nct=nct-1
1008 if(me.eq.king.or..not.out1file)
1009 & write (iout,'(a,i3)') 'nsup=',nsup
1011 if (nsup.le.(nct-nnt+1)) then
1012 do i=0,nct-nnt+1-nsup
1013 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1019 & 'Error - sequences to be superposed do not match.'
1022 do i=0,nsup-(nct-nnt+1)
1023 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1025 nstart_sup=nstart_sup+i
1031 & 'Error - sequences to be superposed do not match.'
1034 if (nsup.eq.0) nsup=nct-nnt
1035 if (nstart_sup.eq.0) nstart_sup=nnt
1036 if (nstart_seq.eq.0) nstart_seq=nnt
1037 if(me.eq.king.or..not.out1file)
1038 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1039 & ' nstart_seq=',nstart_seq
1041 c--- Zscore rms -------
1042 if (nz_start.eq.0) nz_start=nnt
1043 if (nz_end.eq.0 .and. nsup.gt.0) then
1045 else if (nz_end.eq.0) then
1048 if(me.eq.king.or..not.out1file)then
1049 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1050 write (iout,*) 'IZ_SC=',iz_sc
1052 c----------------------
1055 if (.not.pdbref) then
1056 call read_angles(inp,*38)
1058 38 write (iout,'(a)') 'Error reading reference structure.'
1060 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1061 stop 'Error reading reference structure'
1063 39 call chainbuild_extconf
1065 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1072 cref(j,i,kkk)=c(j,i)
1075 call contact(.true.,ncont_ref,icont_ref,co)
1079 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1081 if (constr_dist.gt.0) call read_dist_constr
1082 write (iout,*) "After read_dist_constr nhpb",nhpb
1083 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1085 if(me.eq.king.or..not.out1file)
1086 & write (iout,*) 'Contact order:',co
1088 if(me.eq.king.or..not.out1file)
1089 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1092 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1094 if(me.eq.king.or..not.out1file)
1095 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1096 & icont_ref(1,i),' ',
1097 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1101 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1102 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1103 & modecalc.ne.10) then
1104 C If input structure hasn't been supplied from the PDB file read or generate
1106 if (iranconf.eq.0 .and. .not. extconf) then
1107 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1108 & write (iout,'(a)') 'Initial geometry will be read in.'
1110 read(inp,'(8f10.5)',end=36,err=36)
1111 & ((c(l,k),l=1,3),k=1,nres),
1112 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1113 write (iout,*) "Exit READ_CART"
1114 write (iout,'(8f10.5)')
1115 & ((c(l,k),l=1,3),k=1,nres),
1116 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1117 call int_from_cart1(.true.)
1118 write (iout,*) "Finish INT_TO_CART"
1121 dc(j,i)=c(j,i+1)-c(j,i)
1122 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1126 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1128 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1129 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1135 call read_angles(inp,*36)
1136 call chainbuild_extconf
1139 36 write (iout,'(a)') 'Error reading angle file.'
1141 call mpi_finalize( MPI_COMM_WORLD,IERR )
1143 stop 'Error reading angle file.'
1145 else if (extconf) then
1146 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1147 & write (iout,'(a)') 'Extended chain initial geometry.'
1149 theta(i)=90d0*deg2rad
1152 phi(i)=180d0*deg2rad
1155 alph(i)=110d0*deg2rad
1158 omeg(i)=-120d0*deg2rad
1159 if (itype(i).le.0) omeg(i)=-omeg(i)
1161 call chainbuild_extconf
1163 if(me.eq.king.or..not.out1file)
1164 & write (iout,'(a)') 'Random-generated initial geometry.'
1168 if (me.eq.king .or. fg_rank.eq.0 .and. (
1169 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1173 call gen_rand_conf(itmp,*30)
1175 30 write (iout,*) 'Failed to generate random conformation',
1176 & ', itrial=',itrial
1177 write (*,*) 'Processor:',me,
1178 & ' Failed to generate random conformation',
1188 write (iout,'(a,i3,a)') 'Processor:',me,
1189 & ' error in generating random conformation.'
1190 write (*,'(a,i3,a)') 'Processor:',me,
1191 & ' error in generating random conformation.'
1194 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1199 & ' error in generating random conformation.'
1204 elseif (modecalc.eq.4) then
1205 read (inp,'(a)') intinname
1206 open (intin,file=intinname,status='old',err=333)
1207 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1208 & write (iout,'(a)') 'intinname',intinname
1209 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1211 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1213 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1215 stop 'Error opening angle file.'
1219 C Generate distance constraints, if the PDB structure is to be regularized.
1220 if (nthread.gt.0) then
1221 call read_threadbase
1224 if (me.eq.king .or. .not. out1file)
1226 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1227 write (iout,'(/a,i3,a)')
1228 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1229 write (iout,'(20i4)') (iss(i),i=1,ns)
1231 write(iout,*)"Running with dynamic disulfide-bond formation"
1233 write (iout,'(/a/)') 'Pre-formed links are:'
1239 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1240 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1246 if (ns.gt.0.and.dyn_ss) then
1250 forcon(i-nss)=forcon(i)
1257 dyn_ss_mask(iss(i))=.true.
1260 if (i2ndstr.gt.0) call secstrp2dihc
1261 c call geom_to_var(nvar,x)
1262 c call etotal(energia(0))
1263 c call enerprint(energia(0))
1264 c call briefout(0,etot)
1266 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1267 cd write (iout,'(a)') 'Variable list:'
1268 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1270 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1271 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1272 & 'Processor',myrank,': end reading molecular data.'
1277 c--------------------------------------------------------------------------
1278 logical function seq_comp(itypea,itypeb,length)
1280 integer length,itypea(length),itypeb(length)
1283 if (itypea(i).ne.itypeb(i)) then
1291 c-----------------------------------------------------------------------------
1292 subroutine read_bridge
1293 C Read information about disulfide bridges.
1294 implicit real*8 (a-h,o-z)
1295 include 'DIMENSIONS'
1299 include 'COMMON.IOUNITS'
1300 include 'COMMON.GEO'
1301 include 'COMMON.VAR'
1302 include 'COMMON.INTERACT'
1303 include 'COMMON.LOCAL'
1304 include 'COMMON.NAMES'
1305 include 'COMMON.CHAIN'
1306 include 'COMMON.FFIELD'
1307 include 'COMMON.SBRIDGE'
1308 include 'COMMON.HEADER'
1309 include 'COMMON.CONTROL'
1310 include 'COMMON.DBASE'
1311 include 'COMMON.THREAD'
1312 include 'COMMON.TIME1'
1313 include 'COMMON.SETUP'
1314 C Read bridging residues.
1315 read (inp,*) ns,(iss(i),i=1,ns)
1317 if(me.eq.king.or..not.out1file)
1318 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1319 C Check whether the specified bridging residues are cystines.
1321 if (itype(iss(i)).ne.1) then
1322 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1323 & 'Do you REALLY think that the residue ',
1324 & restyp(itype(iss(i))),i,
1325 & ' can form a disulfide bridge?!!!'
1326 write (*,'(2a,i3,a)')
1327 & 'Do you REALLY think that the residue ',
1328 & restyp(itype(iss(i))),i,
1329 & ' can form a disulfide bridge?!!!'
1331 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1336 C Read preformed bridges.
1338 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1340 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1343 C Check if the residues involved in bridges are in the specified list of
1344 C bridging residues.
1347 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1348 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1349 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1350 & ' contains residues present in other pairs.'
1351 write (*,'(a,i3,a)') 'Disulfide pair',i,
1352 & ' contains residues present in other pairs.'
1354 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1360 if (ihpb(i).eq.iss(j)) goto 10
1362 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1365 if (jhpb(i).eq.iss(j)) goto 20
1367 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1373 ihpb(i)=ihpb(i)+nres
1374 jhpb(i)=jhpb(i)+nres
1380 c----------------------------------------------------------------------------
1381 subroutine read_x(kanal,*)
1382 implicit real*8 (a-h,o-z)
1383 include 'DIMENSIONS'
1384 include 'COMMON.GEO'
1385 include 'COMMON.VAR'
1386 include 'COMMON.CHAIN'
1387 include 'COMMON.IOUNITS'
1388 include 'COMMON.CONTROL'
1389 include 'COMMON.LOCAL'
1390 include 'COMMON.INTERACT'
1391 c Read coordinates from input
1393 read(kanal,'(8f10.5)',end=10,err=10)
1394 & ((c(l,k),l=1,3),k=1,nres),
1395 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1398 c(j,2*nres)=c(j,nres)
1400 call int_from_cart1(.false.)
1403 dc(j,i)=c(j,i+1)-c(j,i)
1404 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1408 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1410 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1411 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1419 c----------------------------------------------------------------------------
1420 subroutine read_threadbase
1421 implicit real*8 (a-h,o-z)
1422 include 'DIMENSIONS'
1423 include 'COMMON.IOUNITS'
1424 include 'COMMON.GEO'
1425 include 'COMMON.VAR'
1426 include 'COMMON.INTERACT'
1427 include 'COMMON.LOCAL'
1428 include 'COMMON.NAMES'
1429 include 'COMMON.CHAIN'
1430 include 'COMMON.FFIELD'
1431 include 'COMMON.SBRIDGE'
1432 include 'COMMON.HEADER'
1433 include 'COMMON.CONTROL'
1434 include 'COMMON.DBASE'
1435 include 'COMMON.THREAD'
1436 include 'COMMON.TIME1'
1437 C Read pattern database for threading.
1438 read (icbase,*) nseq
1440 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1441 & nres_base(2,i),nres_base(3,i)
1442 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1444 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1445 c & nres_base(2,i),nres_base(3,i)
1446 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1450 if (weidis.eq.0.0D0) weidis=0.1D0
1459 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1460 write (iout,'(a,i5)') 'nexcl: ',nexcl
1461 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1464 c------------------------------------------------------------------------------
1465 subroutine setup_var
1466 implicit real*8 (a-h,o-z)
1467 include 'DIMENSIONS'
1468 include 'COMMON.IOUNITS'
1469 include 'COMMON.GEO'
1470 include 'COMMON.VAR'
1471 include 'COMMON.INTERACT'
1472 include 'COMMON.LOCAL'
1473 include 'COMMON.NAMES'
1474 include 'COMMON.CHAIN'
1475 include 'COMMON.FFIELD'
1476 include 'COMMON.SBRIDGE'
1477 include 'COMMON.HEADER'
1478 include 'COMMON.CONTROL'
1479 include 'COMMON.DBASE'
1480 include 'COMMON.THREAD'
1481 include 'COMMON.TIME1'
1482 C Set up variable list.
1488 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1490 ialph(i,1)=nvar+nside
1494 if (indphi.gt.0) then
1496 else if (indback.gt.0) then
1501 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1504 c----------------------------------------------------------------------------
1505 subroutine gen_dist_constr
1506 C Generate CA distance constraints.
1507 implicit real*8 (a-h,o-z)
1508 include 'DIMENSIONS'
1509 include 'COMMON.IOUNITS'
1510 include 'COMMON.GEO'
1511 include 'COMMON.VAR'
1512 include 'COMMON.INTERACT'
1513 include 'COMMON.LOCAL'
1514 include 'COMMON.NAMES'
1515 include 'COMMON.CHAIN'
1516 include 'COMMON.FFIELD'
1517 include 'COMMON.SBRIDGE'
1518 include 'COMMON.HEADER'
1519 include 'COMMON.CONTROL'
1520 include 'COMMON.DBASE'
1521 include 'COMMON.THREAD'
1522 include 'COMMON.TIME1'
1523 dimension itype_pdb(maxres)
1524 common /pizda/ itype_pdb
1526 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1527 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1528 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1530 do i=nstart_sup,nstart_sup+nsup-1
1531 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1532 cd & ' seq_pdb', restyp(itype_pdb(i))
1533 do j=i+2,nstart_sup+nsup-1
1535 ihpb(nhpb)=i+nstart_seq-nstart_sup
1536 jhpb(nhpb)=j+nstart_seq-nstart_sup
1538 dhpb(nhpb)=dist(i,j)
1541 cd write (iout,'(a)') 'Distance constraints:'
1546 cd if (ii.gt.nres) then
1551 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1552 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1553 cd & dhpb(i),forcon(i)
1557 c----------------------------------------------------------------------------
1559 implicit real*8 (a-h,o-z)
1560 include 'DIMENSIONS'
1561 include 'COMMON.MAP'
1562 include 'COMMON.IOUNITS'
1563 character*3 angid(4) /'THE','PHI','ALP','OME'/
1564 character*80 mapcard,ucase
1566 read (inp,'(a)') mapcard
1567 mapcard=ucase(mapcard)
1568 if (index(mapcard,'PHI').gt.0) then
1570 else if (index(mapcard,'THE').gt.0) then
1572 else if (index(mapcard,'ALP').gt.0) then
1574 else if (index(mapcard,'OME').gt.0) then
1577 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1578 stop 'Error - illegal variable spec in MAP card.'
1580 call readi (mapcard,'RES1',res1(imap),0)
1581 call readi (mapcard,'RES2',res2(imap),0)
1582 if (res1(imap).eq.0) then
1583 res1(imap)=res2(imap)
1584 else if (res2(imap).eq.0) then
1585 res2(imap)=res1(imap)
1587 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1589 & 'Error - illegal definition of variable group in MAP.'
1590 stop 'Error - illegal definition of variable group in MAP.'
1592 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1593 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1594 call readi(mapcard,'NSTEP',nstep(imap),0)
1595 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1597 & 'Illegal boundary and/or step size specification in MAP.'
1598 stop 'Illegal boundary and/or step size specification in MAP.'
1603 c----------------------------------------------------------------------------
1605 implicit real*8 (a-h,o-z)
1606 include 'DIMENSIONS'
1607 include 'COMMON.IOUNITS'
1608 include 'COMMON.GEO'
1609 include 'COMMON.CSA'
1610 include 'COMMON.BANK'
1611 include 'COMMON.CONTROL'
1613 character*620 mcmcard
1614 call card_concat(mcmcard)
1616 call readi(mcmcard,'NCONF',nconf,50)
1617 call readi(mcmcard,'NADD',nadd,0)
1618 call readi(mcmcard,'JSTART',jstart,1)
1619 call readi(mcmcard,'JEND',jend,1)
1620 call readi(mcmcard,'NSTMAX',nstmax,500000)
1621 call readi(mcmcard,'N0',n0,1)
1622 call readi(mcmcard,'N1',n1,6)
1623 call readi(mcmcard,'N2',n2,4)
1624 call readi(mcmcard,'N3',n3,0)
1625 call readi(mcmcard,'N4',n4,0)
1626 call readi(mcmcard,'N5',n5,0)
1627 call readi(mcmcard,'N6',n6,10)
1628 call readi(mcmcard,'N7',n7,0)
1629 call readi(mcmcard,'N8',n8,0)
1630 call readi(mcmcard,'N9',n9,0)
1631 call readi(mcmcard,'N14',n14,0)
1632 call readi(mcmcard,'N15',n15,0)
1633 call readi(mcmcard,'N16',n16,0)
1634 call readi(mcmcard,'N17',n17,0)
1635 call readi(mcmcard,'N18',n18,0)
1637 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1639 call readi(mcmcard,'NDIFF',ndiff,2)
1640 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1641 call readi(mcmcard,'IS1',is1,1)
1642 call readi(mcmcard,'IS2',is2,8)
1643 call readi(mcmcard,'NRAN0',nran0,4)
1644 call readi(mcmcard,'NRAN1',nran1,2)
1645 call readi(mcmcard,'IRR',irr,1)
1646 call readi(mcmcard,'NSEED',nseed,20)
1647 call readi(mcmcard,'NTOTAL',ntotal,10000)
1648 call reada(mcmcard,'CUT1',cut1,2.0d0)
1649 call reada(mcmcard,'CUT2',cut2,5.0d0)
1650 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1651 call readi(mcmcard,'ICMAX',icmax,3)
1652 call readi(mcmcard,'IRESTART',irestart,0)
1653 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1656 call reada(mcmcard,'DELE',dele,20.0d0)
1657 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1658 call readi(mcmcard,'IREF',iref,0)
1659 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1660 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1661 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1662 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1663 write (iout,*) "NCONF_IN",nconf_in
1666 c----------------------------------------------------------------------------
1667 cfmc subroutine mcmfread
1668 cfmc implicit real*8 (a-h,o-z)
1669 cfmc include 'DIMENSIONS'
1670 cfmc include 'COMMON.MCMF'
1671 cfmc include 'COMMON.IOUNITS'
1672 cfmc include 'COMMON.GEO'
1673 cfmc character*80 ucase
1674 cfmc character*620 mcmcard
1675 cfmc call card_concat(mcmcard)
1677 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1678 cfmc write(iout,*)'MAXRANT=',maxrant
1679 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1680 cfmc write(iout,*)'MAXFAM=',maxfam
1681 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1682 cfmc write(iout,*)'NNET1=',nnet1
1683 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1684 cfmc write(iout,*)'NNET2=',nnet2
1685 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1686 cfmc write(iout,*)'NNET3=',nnet3
1687 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1688 cfmc write(iout,*)'ILASTT=',ilastt
1689 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1690 cfmc write(iout,*)'MAXSTR=',maxstr
1691 cfmc maxstr_f=maxstr/maxfam
1692 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1693 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1694 cfmc write(iout,*)'NMCMF=',nmcmf
1695 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1696 cfmc write(iout,*)'IFOCUS=',ifocus
1697 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1698 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1699 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1700 cfmc write(iout,*)'INTPRT=',intprt
1701 cfmc call readi(mcmcard,'IPRT',iprt,100)
1702 cfmc write(iout,*)'IPRT=',iprt
1703 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1704 cfmc write(iout,*)'IMAXTR=',imaxtr
1705 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1706 cfmc write(iout,*)'MAXEVEN=',maxeven
1707 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1708 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1709 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1710 cfmc write(iout,*)'INIMIN=',inimin
1711 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1712 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1713 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1714 cfmc write(iout,*)'NTHREAD=',nthread
1715 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1716 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1717 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1718 cfmc write(iout,*)'MAXPERT=',maxpert
1719 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1720 cfmc write(iout,*)'IRMSD=',irmsd
1721 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1722 cfmc write(iout,*)'DENEMIN=',denemin
1723 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1724 cfmc write(iout,*)'RCUT1S=',rcut1s
1725 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1726 cfmc write(iout,*)'RCUT1E=',rcut1e
1727 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1728 cfmc write(iout,*)'RCUT2S=',rcut2s
1729 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1730 cfmc write(iout,*)'RCUT2E=',rcut2e
1731 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1732 cfmc write(iout,*)'DPERT1=',d_pert1
1733 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1734 cfmc write(iout,*)'DPERT1A=',d_pert1a
1735 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1736 cfmc write(iout,*)'DPERT2=',d_pert2
1737 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1738 cfmc write(iout,*)'DPERT2A=',d_pert2a
1739 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1740 cfmc write(iout,*)'DPERT2B=',d_pert2b
1741 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1742 cfmc write(iout,*)'DPERT2C=',d_pert2c
1743 cfmc d_pert1=deg2rad*d_pert1
1744 cfmc d_pert1a=deg2rad*d_pert1a
1745 cfmc d_pert2=deg2rad*d_pert2
1746 cfmc d_pert2a=deg2rad*d_pert2a
1747 cfmc d_pert2b=deg2rad*d_pert2b
1748 cfmc d_pert2c=deg2rad*d_pert2c
1749 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1750 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1751 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1752 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1753 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1754 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1755 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1756 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1757 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1758 cfmc write(iout,*)'RCUTINI=',rcutini
1759 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1760 cfmc write(iout,*)'GRAT=',grat
1761 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1762 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1766 c----------------------------------------------------------------------------
1768 implicit real*8 (a-h,o-z)
1769 include 'DIMENSIONS'
1770 include 'COMMON.MCM'
1771 include 'COMMON.MCE'
1772 include 'COMMON.IOUNITS'
1774 character*320 mcmcard
1775 call card_concat(mcmcard)
1776 call readi(mcmcard,'MAXACC',maxacc,100)
1777 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1778 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1779 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1780 call readi(mcmcard,'MAXREPM',maxrepm,200)
1781 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1782 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1783 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1784 call reada(mcmcard,'E_UP',e_up,5.0D0)
1785 call reada(mcmcard,'DELTE',delte,0.1D0)
1786 call readi(mcmcard,'NSWEEP',nsweep,5)
1787 call readi(mcmcard,'NSTEPH',nsteph,0)
1788 call readi(mcmcard,'NSTEPC',nstepc,0)
1789 call reada(mcmcard,'TMIN',tmin,298.0D0)
1790 call reada(mcmcard,'TMAX',tmax,298.0D0)
1791 call readi(mcmcard,'NWINDOW',nwindow,0)
1792 call readi(mcmcard,'PRINT_MC',print_mc,0)
1793 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1794 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1795 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1796 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1797 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1798 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1799 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1800 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1801 if (nwindow.gt.0) then
1802 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1804 winlen(i)=winend(i)-winstart(i)+1
1807 if (tmax.lt.tmin) tmax=tmin
1808 if (tmax.eq.tmin) then
1812 if (nstepc.gt.0 .and. nsteph.gt.0) then
1813 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1814 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1816 C Probabilities of different move types
1817 sumpro_type(0)=0.0D0
1818 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1819 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1820 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1821 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1822 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1823 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1824 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1826 print *,'i',i,' sumprotype',sumpro_type(i)
1827 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1828 print *,'i',i,' sumprotype',sumpro_type(i)
1832 c----------------------------------------------------------------------------
1833 subroutine read_minim
1834 implicit real*8 (a-h,o-z)
1835 include 'DIMENSIONS'
1836 include 'COMMON.MINIM'
1837 include 'COMMON.IOUNITS'
1839 character*320 minimcard
1840 call card_concat(minimcard)
1841 call readi(minimcard,'MAXMIN',maxmin,2000)
1842 call readi(minimcard,'MAXFUN',maxfun,5000)
1843 call readi(minimcard,'MINMIN',minmin,maxmin)
1844 call readi(minimcard,'MINFUN',minfun,maxmin)
1845 call reada(minimcard,'TOLF',tolf,1.0D-2)
1846 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1847 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1848 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1849 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1850 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1851 & 'Options in energy minimization:'
1852 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1853 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1854 & 'MinMin:',MinMin,' MinFun:',MinFun,
1855 & ' TolF:',TolF,' RTolF:',RTolF
1858 c----------------------------------------------------------------------------
1859 subroutine read_angles(kanal,*)
1860 implicit real*8 (a-h,o-z)
1861 include 'DIMENSIONS'
1862 include 'COMMON.GEO'
1863 include 'COMMON.VAR'
1864 include 'COMMON.CHAIN'
1865 include 'COMMON.IOUNITS'
1866 include 'COMMON.CONTROL'
1867 c Read angles from input
1869 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1870 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1871 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1872 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1875 c 9/7/01 avoid 180 deg valence angle
1876 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1878 theta(i)=deg2rad*theta(i)
1879 phi(i)=deg2rad*phi(i)
1880 alph(i)=deg2rad*alph(i)
1881 omeg(i)=deg2rad*omeg(i)
1886 c----------------------------------------------------------------------------
1887 subroutine reada(rekord,lancuch,wartosc,default)
1889 character*(*) rekord,lancuch
1890 double precision wartosc,default
1893 iread=index(rekord,lancuch)
1894 if (iread.eq.0) then
1898 iread=iread+ilen(lancuch)+1
1899 read (rekord(iread:),*,err=10,end=10) wartosc
1904 c----------------------------------------------------------------------------
1905 subroutine readi(rekord,lancuch,wartosc,default)
1907 character*(*) rekord,lancuch
1908 integer wartosc,default
1911 iread=index(rekord,lancuch)
1912 if (iread.eq.0) then
1916 iread=iread+ilen(lancuch)+1
1917 read (rekord(iread:),*,err=10,end=10) wartosc
1922 c----------------------------------------------------------------------------
1923 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1926 integer tablica(dim),default
1927 character*(*) rekord,lancuch
1934 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1935 if (iread.eq.0) return
1936 iread=iread+ilen(lancuch)+1
1937 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1940 c----------------------------------------------------------------------------
1941 subroutine multreada(rekord,lancuch,tablica,dim,default)
1944 double precision tablica(dim),default
1945 character*(*) rekord,lancuch
1952 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1953 if (iread.eq.0) return
1954 iread=iread+ilen(lancuch)+1
1955 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1958 c----------------------------------------------------------------------------
1959 subroutine openunits
1960 implicit real*8 (a-h,o-z)
1961 include 'DIMENSIONS'
1964 character*16 form,nodename
1967 include 'COMMON.SETUP'
1968 include 'COMMON.IOUNITS'
1970 include 'COMMON.CONTROL'
1971 integer lenpre,lenpot,ilen,lentmp
1973 character*3 out1file_text,ucase
1976 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1977 call getenv_loc("PREFIX",prefix)
1979 call getenv_loc("POT",pot)
1980 call getenv_loc("DIRTMP",tmpdir)
1981 call getenv_loc("CURDIR",curdir)
1982 call getenv_loc("OUT1FILE",out1file_text)
1983 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1984 out1file_text=ucase(out1file_text)
1985 if (out1file_text(1:1).eq."Y") then
1988 out1file=fg_rank.gt.0
1993 if (lentmp.gt.0) then
1994 write (*,'(80(1h!))')
1995 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1996 write (*,'(80(1h!))')
1997 write (*,*)"All output files will be on node /tmp directory."
1999 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2000 if (me.eq.king) then
2001 write (*,*) "The master node is ",nodename
2002 else if (fg_rank.eq.0) then
2003 write (*,*) "I am the CG slave node ",nodename
2005 write (*,*) "I am the FG slave node ",nodename
2008 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2009 lenpre = lentmp+lenpre+1
2011 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2012 C Get the names and open the input files
2013 #if defined(WINIFL) || defined(WINPGI)
2014 open(1,file=pref_orig(:ilen(pref_orig))//
2015 & '.inp',status='old',readonly,shared)
2016 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2017 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2018 C Get parameter filenames and open the parameter files.
2019 call getenv_loc('BONDPAR',bondname)
2020 open (ibond,file=bondname,status='old',readonly,shared)
2021 call getenv_loc('THETPAR',thetname)
2022 open (ithep,file=thetname,status='old',readonly,shared)
2023 call getenv_loc('ROTPAR',rotname)
2024 open (irotam,file=rotname,status='old',readonly,shared)
2025 call getenv_loc('TORPAR',torname)
2026 open (itorp,file=torname,status='old',readonly,shared)
2027 call getenv_loc('TORDPAR',tordname)
2028 open (itordp,file=tordname,status='old',readonly,shared)
2029 call getenv_loc('TORKCC',torkccname)
2030 open (itorkcc,file=torkccname,status='old',readonly,shared)
2031 call getenv_loc('THETKCC',thetkccname)
2032 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2033 call getenv_loc('FOURIER',fouriername)
2034 open (ifourier,file=fouriername,status='old',readonly,shared)
2035 call getenv_loc('ELEPAR',elename)
2036 open (ielep,file=elename,status='old',readonly,shared)
2037 call getenv_loc('SIDEPAR',sidename)
2038 open (isidep,file=sidename,status='old',readonly,shared)
2039 call getenv_loc('LIPTRANPAR',liptranname)
2040 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2041 #elif (defined CRAY) || (defined AIX)
2042 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2044 c print *,"Processor",myrank," opened file 1"
2045 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2046 c print *,"Processor",myrank," opened file 9"
2047 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2048 C Get parameter filenames and open the parameter files.
2049 call getenv_loc('BONDPAR',bondname)
2050 open (ibond,file=bondname,status='old',action='read')
2051 c print *,"Processor",myrank," opened file IBOND"
2052 call getenv_loc('THETPAR',thetname)
2053 open (ithep,file=thetname,status='old',action='read')
2054 c print *,"Processor",myrank," opened file ITHEP"
2055 call getenv_loc('ROTPAR',rotname)
2056 open (irotam,file=rotname,status='old',action='read')
2057 c print *,"Processor",myrank," opened file IROTAM"
2058 call getenv_loc('TORPAR',torname)
2059 open (itorp,file=torname,status='old',action='read')
2060 c print *,"Processor",myrank," opened file ITORP"
2061 call getenv_loc('TORDPAR',tordname)
2062 open (itordp,file=tordname,status='old',action='read')
2063 call getenv_loc('TORKCC',torkccname)
2064 open (itorkcc,file=torkccname,status='old',action='read')
2065 call getenv_loc('THETKCC',thetkccname)
2066 open (ithetkcc,file=thetkccname,status='old',action='read')
2067 c print *,"Processor",myrank," opened file ITORDP"
2068 call getenv_loc('SCCORPAR',sccorname)
2069 open (isccor,file=sccorname,status='old',action='read')
2070 c print *,"Processor",myrank," opened file ISCCOR"
2071 call getenv_loc('FOURIER',fouriername)
2072 open (ifourier,file=fouriername,status='old',action='read')
2073 c print *,"Processor",myrank," opened file IFOURIER"
2074 call getenv_loc('ELEPAR',elename)
2075 open (ielep,file=elename,status='old',action='read')
2076 c print *,"Processor",myrank," opened file IELEP"
2077 call getenv_loc('SIDEPAR',sidename)
2078 open (isidep,file=sidename,status='old',action='read')
2079 call getenv_loc('LIPTRANPAR',liptranname)
2080 open (iliptranpar,file=liptranname,status='old',action='read')
2081 c print *,"Processor",myrank," opened file ISIDEP"
2082 c print *,"Processor",myrank," opened parameter files"
2084 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2085 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2086 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2087 C Get parameter filenames and open the parameter files.
2088 call getenv_loc('BONDPAR',bondname)
2089 open (ibond,file=bondname,status='old')
2090 call getenv_loc('THETPAR',thetname)
2091 open (ithep,file=thetname,status='old')
2092 call getenv_loc('ROTPAR',rotname)
2093 open (irotam,file=rotname,status='old')
2094 call getenv_loc('TORPAR',torname)
2095 open (itorp,file=torname,status='old')
2096 call getenv_loc('TORDPAR',tordname)
2097 open (itordp,file=tordname,status='old')
2098 call getenv_loc('TORKCC',torkccname)
2099 open (itorkcc,file=torkccname,status='old')
2100 call getenv_loc('THETKCC',thetkccname)
2101 open (ithetkcc,file=thetkccname,status='old')
2102 call getenv_loc('SCCORPAR',sccorname)
2103 open (isccor,file=sccorname,status='old')
2104 call getenv_loc('FOURIER',fouriername)
2105 open (ifourier,file=fouriername,status='old')
2106 call getenv_loc('ELEPAR',elename)
2107 open (ielep,file=elename,status='old')
2108 call getenv_loc('SIDEPAR',sidename)
2109 open (isidep,file=sidename,status='old')
2110 call getenv_loc('LIPTRANPAR',liptranname)
2111 open (iliptranpar,file=liptranname,status='old')
2113 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2115 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2116 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2117 C Get parameter filenames and open the parameter files.
2118 call getenv_loc('BONDPAR',bondname)
2119 open (ibond,file=bondname,status='old',readonly)
2120 call getenv_loc('THETPAR',thetname)
2121 open (ithep,file=thetname,status='old',readonly)
2122 call getenv_loc('ROTPAR',rotname)
2123 open (irotam,file=rotname,status='old',readonly)
2124 call getenv_loc('TORPAR',torname)
2125 open (itorp,file=torname,status='old',readonly)
2126 call getenv_loc('TORDPAR',tordname)
2127 open (itordp,file=tordname,status='old',readonly)
2128 call getenv_loc('TORKCC',torkccname)
2129 open (itorkcc,file=torkccname,status='old',readonly)
2130 call getenv_loc('THETKCC',thetkccname)
2131 open (ithetkcc,file=thetkccname,status='old',readonly)
2132 call getenv_loc('SCCORPAR',sccorname)
2133 open (isccor,file=sccorname,status='old',readonly)
2135 call getenv_loc('THETPARPDB',thetname_pdb)
2136 print *,"thetname_pdb ",thetname_pdb
2137 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2138 print *,ithep_pdb," opened"
2140 call getenv_loc('FOURIER',fouriername)
2141 open (ifourier,file=fouriername,status='old',readonly)
2142 call getenv_loc('ELEPAR',elename)
2143 open (ielep,file=elename,status='old',readonly)
2144 call getenv_loc('SIDEPAR',sidename)
2145 open (isidep,file=sidename,status='old',readonly)
2146 call getenv_loc('LIPTRANPAR',liptranname)
2147 open (iliptranpar,file=liptranname,status='old',action='read')
2149 call getenv_loc('ROTPARPDB',rotname_pdb)
2150 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2153 call getenv_loc('TUBEPAR',tubename)
2154 open (itube,file=tubename,status='old',readonly)
2157 C 8/9/01 In the newest version SCp interaction constants are read from a file
2158 C Use -DOLDSCP to use hard-coded constants instead.
2160 call getenv_loc('SCPPAR',scpname)
2161 #if defined(WINIFL) || defined(WINPGI)
2162 open (iscpp,file=scpname,status='old',readonly,shared)
2163 #elif (defined CRAY) || (defined AIX)
2164 open (iscpp,file=scpname,status='old',action='read')
2166 open (iscpp,file=scpname,status='old')
2168 open (iscpp,file=scpname,status='old',readonly)
2171 call getenv_loc('PATTERN',patname)
2172 #if defined(WINIFL) || defined(WINPGI)
2173 open (icbase,file=patname,status='old',readonly,shared)
2174 #elif (defined CRAY) || (defined AIX)
2175 open (icbase,file=patname,status='old',action='read')
2177 open (icbase,file=patname,status='old')
2179 open (icbase,file=patname,status='old',readonly)
2182 C Open output file only for CG processes
2183 c print *,"Processor",myrank," fg_rank",fg_rank
2184 if (fg_rank.eq.0) then
2186 if (nodes.eq.1) then
2189 npos = dlog10(dfloat(nodes-1))+1
2191 if (npos.lt.3) npos=3
2192 write (liczba,'(i1)') npos
2193 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2195 write (liczba,form) me
2196 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2197 & liczba(:ilen(liczba))
2198 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2200 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2202 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2203 & liczba(:ilen(liczba))//'.mol2'
2204 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2205 & liczba(:ilen(liczba))//'.stat'
2207 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2208 & //liczba(:ilen(liczba))//'.stat')
2209 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2212 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2213 & liczba(:ilen(liczba))//'.const'
2218 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2219 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2220 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2221 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2222 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2224 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2226 rest2name=prefix(:ilen(prefix))//'.rst'
2228 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2231 #if defined(AIX) || defined(PGI)
2232 if (me.eq.king .or. .not. out1file)
2233 & open(iout,file=outname,status='unknown')
2235 if (fg_rank.gt.0) then
2236 write (liczba,'(i3.3)') myrank/nfgtasks
2237 write (ll,'(bz,i3.3)') fg_rank
2238 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2243 open(igeom,file=intname,status='unknown',position='append')
2244 open(ipdb,file=pdbname,status='unknown')
2245 open(imol2,file=mol2name,status='unknown')
2246 open(istat,file=statname,status='unknown',position='append')
2248 c1out open(iout,file=outname,status='unknown')
2251 if (me.eq.king .or. .not.out1file)
2252 & open(iout,file=outname,status='unknown')
2254 if (fg_rank.gt.0) then
2255 write (liczba,'(i3.3)') myrank/nfgtasks
2256 write (ll,'(bz,i3.3)') fg_rank
2257 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2262 open(igeom,file=intname,status='unknown',access='append')
2263 open(ipdb,file=pdbname,status='unknown')
2264 open(imol2,file=mol2name,status='unknown')
2265 open(istat,file=statname,status='unknown',access='append')
2267 c1out open(iout,file=outname,status='unknown')
2270 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2271 csa_seed=prefix(:lenpre)//'.CSA.seed'
2272 csa_history=prefix(:lenpre)//'.CSA.history'
2273 csa_bank=prefix(:lenpre)//'.CSA.bank'
2274 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2275 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2276 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2277 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2278 csa_int=prefix(:lenpre)//'.int'
2279 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2280 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2281 csa_in=prefix(:lenpre)//'.CSA.in'
2282 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2285 write (iout,'(80(1h-))')
2286 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2287 write (iout,'(80(1h-))')
2288 write (iout,*) "Input file : ",
2289 & pref_orig(:ilen(pref_orig))//'.inp'
2290 write (iout,*) "Output file : ",
2291 & outname(:ilen(outname))
2293 write (iout,*) "Sidechain potential file : ",
2294 & sidename(:ilen(sidename))
2296 write (iout,*) "SCp potential file : ",
2297 & scpname(:ilen(scpname))
2299 write (iout,*) "Electrostatic potential file : ",
2300 & elename(:ilen(elename))
2301 write (iout,*) "Cumulant coefficient file : ",
2302 & fouriername(:ilen(fouriername))
2303 write (iout,*) "Torsional parameter file : ",
2304 & torname(:ilen(torname))
2305 write (iout,*) "Double torsional parameter file : ",
2306 & tordname(:ilen(tordname))
2307 write (iout,*) "SCCOR parameter file : ",
2308 & sccorname(:ilen(sccorname))
2309 write (iout,*) "Bond & inertia constant file : ",
2310 & bondname(:ilen(bondname))
2311 write (iout,*) "Bending parameter file : ",
2312 & thetname(:ilen(thetname))
2313 write (iout,*) "Rotamer parameter file : ",
2314 & rotname(:ilen(rotname))
2315 write (iout,*) "Thetpdb parameter file : ",
2316 & thetname_pdb(:ilen(thetname_pdb))
2317 write (iout,*) "Threading database : ",
2318 & patname(:ilen(patname))
2320 &write (iout,*)" DIRTMP : ",
2322 write (iout,'(80(1h-))')
2326 c----------------------------------------------------------------------------
2327 subroutine card_concat(card)
2328 implicit real*8 (a-h,o-z)
2329 include 'DIMENSIONS'
2330 include 'COMMON.IOUNITS'
2332 character*80 karta,ucase
2334 read (inp,'(a)') karta
2337 do while (karta(80:80).eq.'&')
2338 card=card(:ilen(card)+1)//karta(:79)
2339 read (inp,'(a)') karta
2342 card=card(:ilen(card)+1)//karta
2345 c----------------------------------------------------------------------------------
2347 implicit real*8 (a-h,o-z)
2348 include 'DIMENSIONS'
2349 include 'COMMON.CHAIN'
2350 include 'COMMON.IOUNITS'
2352 open(irest2,file=rest2name,status='unknown')
2353 read(irest2,*) totT,EK,potE,totE,t_bath
2356 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2359 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2362 read (irest2,*) iset
2367 c---------------------------------------------------------------------------------
2368 subroutine read_fragments
2369 implicit real*8 (a-h,o-z)
2370 include 'DIMENSIONS'
2374 include 'COMMON.SETUP'
2375 include 'COMMON.CHAIN'
2376 include 'COMMON.IOUNITS'
2378 include 'COMMON.CONTROL'
2379 read(inp,*) nset,nfrag,npair,nfrag_back
2380 if(me.eq.king.or..not.out1file)
2381 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2382 & " nfrag_back",nfrag_back
2384 read(inp,*) mset(iset)
2386 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2388 if(me.eq.king.or..not.out1file)
2389 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2390 & ifrag(2,i,iset), qinfrag(i,iset)
2393 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2395 if(me.eq.king.or..not.out1file)
2396 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2397 & ipair(2,i,iset), qinpair(i,iset)
2400 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2401 & wfrag_back(3,i,iset),
2402 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2403 if(me.eq.king.or..not.out1file)
2404 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2405 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2410 C---------------------------------------------------------------------------
2411 subroutine read_afminp
2412 implicit real*8 (a-h,o-z)
2413 include 'DIMENSIONS'
2417 include 'COMMON.SETUP'
2418 include 'COMMON.CONTROL'
2419 include 'COMMON.CHAIN'
2420 include 'COMMON.IOUNITS'
2421 include 'COMMON.SBRIDGE'
2422 character*320 afmcard
2424 call card_concat(afmcard)
2425 call readi(afmcard,"BEG",afmbeg,0)
2426 call readi(afmcard,"END",afmend,0)
2427 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2428 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2429 print *,'FORCE=' ,forceAFMconst
2430 CCCC NOW PROPERTIES FOR AFM
2433 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2435 distafminit=dsqrt(distafminit)
2436 print *,'initdist',distafminit
2440 c-------------------------------------------------------------------------------
2441 subroutine read_dist_constr
2442 implicit real*8 (a-h,o-z)
2443 include 'DIMENSIONS'
2447 include 'COMMON.SETUP'
2448 include 'COMMON.CONTROL'
2449 include 'COMMON.CHAIN'
2450 include 'COMMON.IOUNITS'
2451 include 'COMMON.SBRIDGE'
2452 integer ifrag_(2,100),ipair_(2,100)
2453 double precision wfrag_(100),wpair_(100)
2454 character*500 controlcard
2456 write (iout,*) "Calling read_dist_constr"
2457 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2459 call card_concat(controlcard)
2460 call readi(controlcard,"NFRAG",nfrag_,0)
2461 call readi(controlcard,"NPAIR",npair_,0)
2462 call readi(controlcard,"NDIST",ndist_,0)
2463 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2464 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2465 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2466 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2467 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2468 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2469 c write (iout,*) "IFRAG"
2471 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2473 c write (iout,*) "IPAIR"
2475 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2479 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2480 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2481 & ifrag_(2,i)=nstart_sup+nsup-1
2482 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2484 if (wfrag_(i).gt.0.0d0) then
2485 do j=ifrag_(1,i),ifrag_(2,i)-1
2486 do k=j+1,ifrag_(2,i)
2487 c write (iout,*) "j",j," k",k
2489 if (constr_dist.eq.1) then
2494 forcon(nhpb)=wfrag_(i)
2495 else if (constr_dist.eq.2) then
2496 if (ddjk.le.dist_cut) then
2501 forcon(nhpb)=wfrag_(i)
2508 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2511 if (.not.out1file .or. me.eq.king)
2512 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2513 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2515 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2516 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2523 if (wpair_(i).gt.0.0d0) then
2531 do j=ifrag_(1,ii),ifrag_(2,ii)
2532 do k=ifrag_(1,jj),ifrag_(2,jj)
2536 forcon(nhpb)=wpair_(i)
2537 dhpb(nhpb)=dist(j,k)
2539 if (.not.out1file .or. me.eq.king)
2540 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2541 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2543 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2544 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2552 if (constr_dist.eq.11) then
2553 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2554 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2555 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2558 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2559 & ibecarb(i),forcon(nhpb+1)
2561 if (forcon(nhpb+1).gt.0.0d0) then
2563 if (ibecarb(i).gt.0) then
2564 ihpb(i)=ihpb(i)+nres
2565 jhpb(i)=jhpb(i)+nres
2567 if (dhpb(nhpb).eq.0.0d0)
2568 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2570 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2571 C if (forcon(nhpb+1).gt.0.0d0) then
2573 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2575 if (.not.out1file .or. me.eq.king)
2576 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2577 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2579 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2580 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2587 c-------------------------------------------------------------------------------
2589 subroutine flush(iu)
2594 subroutine flush(iu)
2599 c------------------------------------------------------------------------------
2600 subroutine copy_to_tmp(source)
2601 include "DIMENSIONS"
2602 include "COMMON.IOUNITS"
2603 character*(*) source
2604 character* 256 tmpfile
2608 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2609 inquire(file=tmpfile,exist=ex)
2611 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2612 & " to temporary directory..."
2613 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2614 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2618 c------------------------------------------------------------------------------
2619 subroutine move_from_tmp(source)
2620 include "DIMENSIONS"
2621 include "COMMON.IOUNITS"
2622 character*(*) source
2625 write (*,*) "Moving ",source(:ilen(source)),
2626 & " from temporary directory to working directory"
2627 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2628 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2631 c------------------------------------------------------------------------------
2632 subroutine random_init(seed)
2634 C Initialize random number generator
2636 implicit real*8 (a-h,o-z)
2637 include 'DIMENSIONS'
2640 logical OKRandom, prng_restart
2642 integer iseed_array(4)
2644 include 'COMMON.IOUNITS'
2645 include 'COMMON.TIME1'
2646 include 'COMMON.THREAD'
2647 include 'COMMON.SBRIDGE'
2648 include 'COMMON.CONTROL'
2649 include 'COMMON.MCM'
2650 include 'COMMON.MAP'
2651 include 'COMMON.HEADER'
2652 include 'COMMON.CSA'
2653 include 'COMMON.CHAIN'
2654 include 'COMMON.MUCA'
2656 include 'COMMON.FFIELD'
2657 include 'COMMON.SETUP'
2658 iseed=-dint(dabs(seed))
2659 if (iseed.eq.0) then
2660 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2661 & 'Random seed undefined. The program will stop.'
2662 write (*,'(/80(1h*)/20x,a/80(1h*))')
2663 & 'Random seed undefined. The program will stop.'
2665 call mpi_finalize(mpi_comm_world,ierr)
2667 stop 'Bad random seed.'
2670 if (fg_rank.eq.0) then
2674 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2675 OKRandom = prng_restart(me,iseed)
2678 tmp=65536.0d0**(4-i)
2679 iseed_array(i) = dint(seed/tmp)
2680 seed=seed-iseed_array(i)*tmp
2683 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2684 & (iseed_array(i),i=1,4)
2685 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2686 & (iseed_array(i),i=1,4)
2687 OKRandom = prng_restart(me,iseed_array)
2690 c r1 = prng_next(me)
2691 r1=ran_number(0.0D0,1.0D0)
2693 & write (iout,*) 'ran_num',r1
2694 if (r1.lt.0.0d0) OKRandom=.false.
2696 if (.not.OKRandom) then
2697 write (iout,*) 'PRNG IS NOT WORKING!!!'
2698 print *,'PRNG IS NOT WORKING!!!'
2701 call mpi_abort(mpi_comm_world,error_msg,ierr)
2704 write (iout,*) 'too many processors for parallel prng'
2705 write (*,*) 'too many processors for parallel prng'
2713 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)