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,'GENCONSTR',genconstr,0)
151 C write (iout,*) TUBElog,"TUBEMODE"
152 call readi(controlcard,'IPRINT',iprint,0)
153 C SHIELD keyword sets if the shielding effect of side-chains is used
154 C 0 denots no shielding is used all peptide are equally despite the
155 C solvent accesible area
156 C 1 the newly introduced function
157 C 2 reseved for further possible developement
158 call readi(controlcard,'SHIELD',shield_mode,0)
159 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
160 write(iout,*) "shield_mode",shield_mode
162 call readi(controlcard,'TORMODE',tor_mode,0)
163 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
164 write(iout,*) "torsional and valence angle mode",tor_mode
165 call readi(controlcard,'MAXGEN',maxgen,10000)
166 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
167 call readi(controlcard,"KDIAG",kdiag,0)
168 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
169 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
170 & write (iout,*) "RESCALE_MODE",rescale_mode
171 split_ene=index(controlcard,'SPLIT_ENE').gt.0
172 if (index(controlcard,'REGULAR').gt.0.0D0) then
173 call reada(controlcard,'WEIDIS',weidis,0.1D0)
177 if (index(controlcard,'CHECKGRAD').gt.0) then
179 if (index(controlcard,'CART').gt.0) then
181 elseif (index(controlcard,'CARINT').gt.0) then
186 call reada(controlcard,'DELTA',aincr,1.0d-7)
187 elseif (index(controlcard,'THREAD').gt.0) then
189 call readi(controlcard,'THREAD',nthread,0)
190 if (nthread.gt.0) then
191 call reada(controlcard,'WEIDIS',weidis,0.1D0)
194 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
195 stop 'Error termination in Read_Control.'
197 else if (index(controlcard,'MCMA').gt.0) then
199 else if (index(controlcard,'MCEE').gt.0) then
201 else if (index(controlcard,'MULTCONF').gt.0) then
203 else if (index(controlcard,'MAP').gt.0) then
205 call readi(controlcard,'MAP',nmap,0)
206 else if (index(controlcard,'CSA').gt.0) then
208 crc else if (index(controlcard,'ZSCORE').gt.0) then
210 crc ZSCORE is rm from UNRES, modecalc=9 is available
213 cfcm else if (index(controlcard,'MCMF').gt.0) then
215 else if (index(controlcard,'SOFTREG').gt.0) then
217 else if (index(controlcard,'CHECK_BOND').gt.0) then
219 else if (index(controlcard,'TEST').gt.0) then
221 else if (index(controlcard,'MD').gt.0) then
223 else if (index(controlcard,'RE ').gt.0) then
227 lmuca=index(controlcard,'MUCA').gt.0
228 call readi(controlcard,'MUCADYN',mucadyn,0)
229 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
230 if (lmuca .and. (me.eq.king .or. .not.out1file ))
232 write (iout,*) 'MUCADYN=',mucadyn
233 write (iout,*) 'MUCASMOOTH=',muca_smooth
236 iscode=index(controlcard,'ONE_LETTER')
237 indphi=index(controlcard,'PHI')
238 indback=index(controlcard,'BACK')
239 iranconf=index(controlcard,'RAND_CONF')
240 i2ndstr=index(controlcard,'USE_SEC_PRED')
241 gradout=index(controlcard,'GRADOUT').gt.0
242 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
243 C DISTCHAINMAX become obsolete for periodic boundry condition
244 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
245 C Reading the dimensions of box in x,y,z coordinates
246 call reada(controlcard,'BOXX',boxxsize,100.0d0)
247 call reada(controlcard,'BOXY',boxysize,100.0d0)
248 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
249 c Cutoff range for interactions
250 call reada(controlcard,"R_CUT",r_cut,15.0d0)
251 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
252 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
253 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
254 if (lipthick.gt.0.0d0) then
255 bordliptop=(boxzsize+lipthick)/2.0
256 bordlipbot=bordliptop-lipthick
258 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
259 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
260 buflipbot=bordlipbot+lipbufthick
261 bufliptop=bordliptop-lipbufthick
262 if ((lipbufthick*2.0d0).gt.lipthick)
263 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
265 write(iout,*) "bordliptop=",bordliptop
266 write(iout,*) "bordlipbot=",bordlipbot
267 write(iout,*) "bufliptop=",bufliptop
268 write(iout,*) "buflipbot=",buflipbot
269 write (iout,*) "SHIELD MODE",shield_mode
270 if (TUBElog.gt.0) then
271 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
272 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
273 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
274 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
275 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
276 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
277 buftubebot=bordtubebot+tubebufthick
278 buftubetop=bordtubetop-tubebufthick
280 if (shield_mode.gt.0) then
282 C VSolvSphere the volume of solving sphere
284 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
285 C there will be no distinction between proline peptide group and normal peptide
286 C group in case of shielding parameters
287 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
288 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
289 write (iout,*) VSolvSphere,VSolvSphere_div
290 C long axis of side chain
292 long_r_sidechain(i)=vbldsc0(1,i)
293 short_r_sidechain(i)=sigma0(i)
297 if (me.eq.king .or. .not.out1file )
298 & write (iout,*) "DISTCHAINMAX",distchainmax
300 if(me.eq.king.or..not.out1file)
301 & write (iout,'(2a)') diagmeth(kdiag),
302 & ' routine used to diagonalize matrices.'
305 c--------------------------------------------------------------------------
306 subroutine read_REMDpar
310 implicit real*8 (a-h,o-z)
312 include 'COMMON.IOUNITS'
313 include 'COMMON.TIME1'
316 include 'COMMON.LANGEVIN'
318 include 'COMMON.LANGEVIN.lang0'
320 include 'COMMON.INTERACT'
321 include 'COMMON.NAMES'
323 include 'COMMON.REMD'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SETUP'
327 character*320 controlcard
328 character*3200 controlcard1
329 integer iremd_m_total
331 if(me.eq.king.or..not.out1file)
332 & write (iout,*) "REMD setup"
334 call card_concat(controlcard)
335 call readi(controlcard,"NREP",nrep,3)
336 call readi(controlcard,"NSTEX",nstex,1000)
337 call reada(controlcard,"RETMIN",retmin,10.0d0)
338 call reada(controlcard,"RETMAX",retmax,1000.0d0)
339 mremdsync=(index(controlcard,'SYNC').gt.0)
340 call readi(controlcard,"NSYN",i_sync_step,100)
341 restart1file=(index(controlcard,'REST1FILE').gt.0)
342 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
343 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
344 if(max_cache_traj_use.gt.max_cache_traj)
345 & max_cache_traj_use=max_cache_traj
346 if(me.eq.king.or..not.out1file) then
347 cd if (traj1file) then
348 crc caching is in testing - NTWX is not ignored
349 cd write (iout,*) "NTWX value is ignored"
350 cd write (iout,*) " trajectory is stored to one file by master"
351 cd write (iout,*) " before exchange at NSTEX intervals"
353 write (iout,*) "NREP= ",nrep
354 write (iout,*) "NSTEX= ",nstex
355 write (iout,*) "SYNC= ",mremdsync
356 write (iout,*) "NSYN= ",i_sync_step
357 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
360 if (index(controlcard,'TLIST').gt.0) then
362 call card_concat(controlcard1)
363 read(controlcard1,*) (remd_t(i),i=1,nrep)
364 if(me.eq.king.or..not.out1file)
365 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
368 if (index(controlcard,'MLIST').gt.0) then
370 call card_concat(controlcard1)
371 read(controlcard1,*) (remd_m(i),i=1,nrep)
372 if(me.eq.king.or..not.out1file) then
373 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
376 iremd_m_total=iremd_m_total+remd_m(i)
378 write (iout,*) 'Total number of replicas ',iremd_m_total
381 if(me.eq.king.or..not.out1file)
382 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
385 c--------------------------------------------------------------------------
386 subroutine read_MDpar
390 implicit real*8 (a-h,o-z)
392 include 'COMMON.IOUNITS'
393 include 'COMMON.TIME1'
396 include 'COMMON.LANGEVIN'
398 include 'COMMON.LANGEVIN.lang0'
400 include 'COMMON.INTERACT'
401 include 'COMMON.NAMES'
403 include 'COMMON.SETUP'
404 include 'COMMON.CONTROL'
405 include 'COMMON.SPLITELE'
407 character*320 controlcard
409 call card_concat(controlcard)
410 call readi(controlcard,"NSTEP",n_timestep,1000000)
411 call readi(controlcard,"NTWE",ntwe,100)
412 call readi(controlcard,"NTWX",ntwx,1000)
413 call reada(controlcard,"DT",d_time,1.0d-1)
414 call reada(controlcard,"DVMAX",dvmax,2.0d1)
415 call reada(controlcard,"DAMAX",damax,1.0d1)
416 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
417 call readi(controlcard,"LANG",lang,0)
418 RESPA = index(controlcard,"RESPA") .gt. 0
419 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
420 ntime_split0=ntime_split
421 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
422 ntime_split0=ntime_split
423 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
424 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
425 rest = index(controlcard,"REST").gt.0
426 tbf = index(controlcard,"TBF").gt.0
427 usampl = index(controlcard,"USAMPL").gt.0
429 mdpdb = index(controlcard,"MDPDB").gt.0
430 call reada(controlcard,"T_BATH",t_bath,300.0d0)
431 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
432 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
433 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
434 if (count_reset_moment.eq.0) count_reset_moment=1000000000
435 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
436 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
437 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
438 if (count_reset_vel.eq.0) count_reset_vel=1000000000
439 large = index(controlcard,"LARGE").gt.0
440 print_compon = index(controlcard,"PRINT_COMPON").gt.0
441 rattle = index(controlcard,"RATTLE").gt.0
442 c if performing umbrella sampling, fragments constrained are read from the fragment file
448 if(me.eq.king.or..not.out1file) then
450 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
452 write (iout,'(a)') "The units are:"
453 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
454 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
455 & " acceleration: angstrom/(48.9 fs)**2"
456 write (iout,'(a)') "energy: kcal/mol, temperature: K"
458 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
459 write (iout,'(a60,f10.5,a)')
460 & "Initial time step of numerical integration:",d_time,
462 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
464 write (iout,'(2a,i4,a)')
465 & "A-MTS algorithm used; initial time step for fast-varying",
466 & " short-range forces split into",ntime_split," steps."
467 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
468 & r_cut," lambda",rlamb
470 write (iout,'(2a,f10.5)')
471 & "Maximum acceleration threshold to reduce the time step",
472 & "/increase split number:",damax
473 write (iout,'(2a,f10.5)')
474 & "Maximum predicted energy drift to reduce the timestep",
475 & "/increase split number:",edriftmax
476 write (iout,'(a60,f10.5)')
477 & "Maximum velocity threshold to reduce velocities:",dvmax
478 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
479 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
480 if (rattle) write (iout,'(a60)')
481 & "Rattle algorithm used to constrain the virtual bonds"
485 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
486 call reada(controlcard,"RWAT",rwat,1.4d0)
487 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
488 surfarea=index(controlcard,"SURFAREA").gt.0
489 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
490 if(me.eq.king.or..not.out1file)then
491 write (iout,'(/a,$)') "Langevin dynamics calculation"
494 & " with direct integration of Langevin equations"
495 else if (lang.eq.2) then
496 write (iout,'(a/)') " with TINKER stochasic MD integrator"
497 else if (lang.eq.3) then
498 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
499 else if (lang.eq.4) then
500 write (iout,'(a/)') " in overdamped mode"
502 write (iout,'(//a,i5)')
503 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
506 write (iout,'(a60,f10.5)') "Temperature:",t_bath
507 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
508 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
509 write (iout,'(a60,f10.5)')
510 & "Scaling factor of the friction forces:",scal_fric
511 if (surfarea) write (iout,'(2a,i10,a)')
512 & "Friction coefficients will be scaled by solvent-accessible",
513 & " surface area every",reset_fricmat," steps."
515 c Calculate friction coefficients and bounds of stochastic forces
516 eta=6*pi*cPoise*etawat
517 if(me.eq.king.or..not.out1file)
518 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
520 gamp=scal_fric*(pstok+rwat)*eta
521 stdfp=dsqrt(2*Rb*t_bath/d_time)
523 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
524 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
526 if(me.eq.king.or..not.out1file)then
527 write (iout,'(/2a/)')
528 & "Radii of site types and friction coefficients and std's of",
529 & " stochastic forces of fully exposed sites"
530 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
532 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
533 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
537 if(me.eq.king.or..not.out1file)then
538 write (iout,'(a)') "Berendsen bath calculation"
539 write (iout,'(a60,f10.5)') "Temperature:",t_bath
540 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
542 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
543 & count_reset_moment," steps"
545 & write (iout,'(a,i10,a)')
546 & "Velocities will be reset at random every",count_reset_vel,
550 if(me.eq.king.or..not.out1file)
551 & write (iout,'(a31)') "Microcanonical mode calculation"
553 if(me.eq.king.or..not.out1file)then
554 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
556 write(iout,*) "MD running with constraints."
557 write(iout,*) "Equilibration time ", eq_time, " mtus."
558 write(iout,*) "Constraining ", nfrag," fragments."
559 write(iout,*) "Length of each fragment, weight and q0:"
561 write (iout,*) "Set of restraints #",iset
563 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
564 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
566 write(iout,*) "constraints between ", npair, "fragments."
567 write(iout,*) "constraint pairs, weights and q0:"
569 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
570 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
572 write(iout,*) "angle constraints within ", nfrag_back,
573 & "backbone fragments."
574 write(iout,*) "fragment, weights:"
576 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
577 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
578 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
581 iset=mod(kolor,nset)+1
584 if(me.eq.king.or..not.out1file)
585 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
588 c------------------------------------------------------------------------------
591 C Read molecular data.
593 implicit real*8 (a-h,o-z)
599 include 'COMMON.IOUNITS'
602 include 'COMMON.INTERACT'
603 include 'COMMON.LOCAL'
604 include 'COMMON.NAMES'
605 include 'COMMON.CHAIN'
606 include 'COMMON.FFIELD'
607 include 'COMMON.SBRIDGE'
608 include 'COMMON.HEADER'
609 include 'COMMON.CONTROL'
610 include 'COMMON.DBASE'
611 include 'COMMON.THREAD'
612 include 'COMMON.CONTACTS'
613 include 'COMMON.TORCNSTR'
614 include 'COMMON.TIME1'
615 include 'COMMON.BOUNDS'
617 include 'COMMON.SETUP'
618 include 'COMMON.SHIELD'
619 character*4 sequence(maxres)
621 double precision x(maxvar)
622 character*256 pdbfile
623 character*400 weightcard
624 character*80 weightcard_t,ucase
625 dimension itype_pdb(maxres)
626 common /pizda/ itype_pdb
627 logical seq_comp,fail
628 double precision energia(0:n_ene)
634 C Read weights of the subsequent energy terms.
635 call card_concat(weightcard)
636 call reada(weightcard,'WLONG',wlong,1.0D0)
637 call reada(weightcard,'WSC',wsc,wlong)
638 call reada(weightcard,'WSCP',wscp,wlong)
639 call reada(weightcard,'WELEC',welec,1.0D0)
640 call reada(weightcard,'WVDWPP',wvdwpp,welec)
641 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
642 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
643 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
644 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
645 call reada(weightcard,'WTURN3',wturn3,1.0D0)
646 call reada(weightcard,'WTURN4',wturn4,1.0D0)
647 call reada(weightcard,'WTURN6',wturn6,1.0D0)
648 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
649 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
650 call reada(weightcard,'WBOND',wbond,1.0D0)
651 call reada(weightcard,'WTOR',wtor,1.0D0)
652 call reada(weightcard,'WTORD',wtor_d,1.0D0)
653 call reada(weightcard,'WANG',wang,1.0D0)
654 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
655 call reada(weightcard,'SCAL14',scal14,0.4D0)
656 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
657 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
658 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
659 call reada(weightcard,'TEMP0',temp0,300.0d0)
660 call reada(weightcard,'WSHIELD',wshield,1.0d0)
661 call reada(weightcard,'WLT',wliptran,0.0D0)
662 call reada(weightcard,'WTUBE',wtube,1.0D0)
663 if (index(weightcard,'SOFT').gt.0) ipot=6
664 C 12/1/95 Added weight for the multi-body term WCORR
665 call reada(weightcard,'WCORRH',wcorr,1.0D0)
666 if (wcorr4.gt.0.0d0) wcorr=wcorr4
686 if(me.eq.king.or..not.out1file)
687 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
688 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
690 10 format (/'Energy-term weights (unscaled):'//
691 & 'WSCC= ',f10.6,' (SC-SC)'/
692 & 'WSCP= ',f10.6,' (SC-p)'/
693 & 'WELEC= ',f10.6,' (p-p electr)'/
694 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
695 & 'WBOND= ',f10.6,' (stretching)'/
696 & 'WANG= ',f10.6,' (bending)'/
697 & 'WSCLOC= ',f10.6,' (SC local)'/
698 & 'WTOR= ',f10.6,' (torsional)'/
699 & 'WTORD= ',f10.6,' (double torsional)'/
700 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
701 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
702 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
703 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
704 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
705 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
706 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
707 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
708 & 'WTURN6= ',f10.6,' (turns, 6th order)')
709 if(me.eq.king.or..not.out1file)then
710 if (wcorr4.gt.0.0d0) then
711 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
712 & 'between contact pairs of peptide groups'
713 write (iout,'(2(a,f5.3/))')
714 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
715 & 'Range of quenching the correlation terms:',2*delt_corr
716 else if (wcorr.gt.0.0d0) then
717 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
718 & 'between contact pairs of peptide groups'
720 write (iout,'(a,f8.3)')
721 & 'Scaling factor of 1,4 SC-p interactions:',scal14
722 write (iout,'(a,f8.3)')
723 & 'General scaling factor of SC-p interactions:',scalscp
725 r0_corr=cutoff_corr-delt_corr
727 aad(i,1)=scalscp*aad(i,1)
728 aad(i,2)=scalscp*aad(i,2)
729 bad(i,1)=scalscp*bad(i,1)
730 bad(i,2)=scalscp*bad(i,2)
732 call rescale_weights(t_bath)
733 if(me.eq.king.or..not.out1file)
734 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
735 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
737 22 format (/'Energy-term weights (scaled):'//
738 & 'WSCC= ',f10.6,' (SC-SC)'/
739 & 'WSCP= ',f10.6,' (SC-p)'/
740 & 'WELEC= ',f10.6,' (p-p electr)'/
741 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
742 & 'WBOND= ',f10.6,' (stretching)'/
743 & 'WANG= ',f10.6,' (bending)'/
744 & 'WSCLOC= ',f10.6,' (SC local)'/
745 & 'WTOR= ',f10.6,' (torsional)'/
746 & 'WTORD= ',f10.6,' (double torsional)'/
747 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
748 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
749 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
750 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
751 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
752 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
753 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
754 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
755 & 'WTURN6= ',f10.6,' (turns, 6th order)')
756 if(me.eq.king.or..not.out1file)
757 & write (iout,*) "Reference temperature for weights calculation:",
759 call reada(weightcard,"D0CM",d0cm,3.78d0)
760 call reada(weightcard,"AKCM",akcm,15.1d0)
761 call reada(weightcard,"AKTH",akth,11.0d0)
762 call reada(weightcard,"AKCT",akct,12.0d0)
763 call reada(weightcard,"V1SS",v1ss,-1.08d0)
764 call reada(weightcard,"V2SS",v2ss,7.61d0)
765 call reada(weightcard,"V3SS",v3ss,13.7d0)
766 call reada(weightcard,"EBR",ebr,-5.50D0)
767 call reada(weightcard,"ATRISS",atriss,0.301D0)
768 call reada(weightcard,"BTRISS",btriss,0.021D0)
769 call reada(weightcard,"CTRISS",ctriss,1.001D0)
770 call reada(weightcard,"DTRISS",dtriss,1.001D0)
771 call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
772 write (iout,*) "ATRISS=", atriss
773 write (iout,*) "BTRISS=", btriss
774 write (iout,*) "CTRISS=", ctriss
775 write (iout,*) "DTRISS=", dtriss
776 if (shield_mode.gt.0) then
778 write (iout,*) "liscale not used in shield mode"
780 write (iout,*) "lipid scaling factor", lipscale
781 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
783 dyn_ss_mask(i)=.false.
787 dyn_ssbond_ij(i,j)=1.0d300
790 call reada(weightcard,"HT",Ht,0.0D0)
792 ss_depth=ebr/wsc-0.25*eps(1,1)
793 Ht=Ht/wsc-0.25*eps(1,1)
801 if (wstrain.ne.0.0) then
802 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
807 write (iout,*) "wshield,", wshield
808 if(me.eq.king.or..not.out1file) then
809 write (iout,*) "Parameters of the SS-bond potential:"
810 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
812 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
813 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
814 write (iout,*)" HT",Ht
815 print *,'indpdb=',indpdb,' pdbref=',pdbref
817 if (indpdb.gt.0 .or. pdbref) then
818 read(inp,'(a)') pdbfile
819 if(me.eq.king.or..not.out1file)
820 & write (iout,'(2a)') 'PDB data will be read from file ',
821 & pdbfile(:ilen(pdbfile))
822 open(ipdbin,file=pdbfile,status='old',err=33)
824 33 write (iout,'(a)') 'Error opening PDB file.'
827 c write (iout,*) 'Begin reading pdb data'
830 c write (iout,*) 'Finished reading pdb data'
832 if(me.eq.king.or..not.out1file)
833 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
834 & ' nstart_sup=',nstart_sup
836 itype_pdb(i)=itype(i)
840 nct=nstart_sup+nsup-1
841 call contact(.false.,ncont_ref,icont_ref,co)
844 if(me.eq.king.or..not.out1file)
845 & write(iout,*)'Adding sidechains'
849 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
852 do while (fail.and.nsi.le.maxsi)
853 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
856 if(fail) write(iout,*)'Adding sidechain failed for res ',
857 & i,' after ',nsi,' trials'
862 if (indpdb.eq.0) then
863 C Read sequence if not taken from the pdb file.
865 c print *,'nres=',nres
866 if (iscode.gt.0) then
867 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
869 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
871 C Convert sequence to numeric code
873 itype(i)=rescode(i,sequence(i),iscode)
875 C Assign initial virtual bond lengths
881 vbld(i+nres)=dsc(iabs(itype(i)))
882 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
883 c write (iout,*) "i",i," itype",itype(i),
884 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
888 c print '(20i4)',(itype(i),i=1,nres)
891 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
893 if (itype(i).eq.ntyp1) then
897 else if (iabs(itype(i+1)).ne.20) then
899 else if (iabs(itype(i)).ne.20) then
906 if(me.eq.king.or..not.out1file)then
907 write (iout,*) "ITEL"
909 write (iout,*) i,itype(i),itel(i)
911 print *,'Call Read_Bridge.'
914 C 8/13/98 Set limits to generating the dihedral angles
919 read (inp,*) ndih_constr
920 if (ndih_constr.gt.0) then
922 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
924 if(me.eq.king.or..not.out1file)then
926 & 'There are',ndih_constr,' constraints on phi angles.'
928 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
933 phi0(i)=deg2rad*phi0(i)
934 drange(i)=deg2rad*drange(i)
936 C if(me.eq.king.or..not.out1file)
937 C & write (iout,*) 'FTORS',ftors
940 phibound(1,ii) = phi0(i)-drange(i)
941 phibound(2,ii) = phi0(i)+drange(i)
944 C first setting the theta boundaries to 0 to pi
945 C this mean that there is no energy penalty for any angle occuring this can be applied
946 C for generate random conformation but is not implemented in this way
951 C begin reading theta constrains this is quartic constrains allowing to
952 C have smooth second derivative
953 if (with_theta_constr) then
954 C with_theta_constr is keyword allowing for occurance of theta constrains
955 read (inp,*) ntheta_constr
956 C ntheta_constr is the number of theta constrains
957 if (ntheta_constr.gt.0) then
959 read (inp,*) (itheta_constr(i),theta_constr0(i),
960 & theta_drange(i),for_thet_constr(i),
962 C the above code reads from 1 to ntheta_constr
963 C itheta_constr(i) residue i for which is theta_constr
964 C theta_constr0 the global minimum value
965 C theta_drange is range for which there is no energy penalty
966 C for_thet_constr is the force constant for quartic energy penalty
968 if(me.eq.king.or..not.out1file)then
970 & 'There are',ntheta_constr,' constraints on phi angles.'
972 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
978 theta_constr0(i)=deg2rad*theta_constr0(i)
979 theta_drange(i)=deg2rad*theta_drange(i)
981 C if(me.eq.king.or..not.out1file)
982 C & write (iout,*) 'FTORS',ftors
983 C do i=1,ntheta_constr
984 C ii = itheta_constr(i)
985 C thetabound(1,ii) = phi0(i)-drange(i)
986 C thetabound(2,ii) = phi0(i)+drange(i)
988 endif ! ntheta_constr.gt.0
989 endif! with_theta_constr
991 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
992 C write (iout,*) "with_dihed_constr ",with_dihed_constr
997 write (iout,'(a)') 'Boundaries in phi angle sampling:'
999 write (iout,'(a3,i5,2f10.1)')
1000 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1006 cd print *,'NNT=',NNT,' NCT=',NCT
1007 if (itype(1).eq.ntyp1) nnt=2
1008 if (itype(nres).eq.ntyp1) nct=nct-1
1010 if(me.eq.king.or..not.out1file)
1011 & write (iout,'(a,i3)') 'nsup=',nsup
1013 if (nsup.le.(nct-nnt+1)) then
1014 do i=0,nct-nnt+1-nsup
1015 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1021 & 'Error - sequences to be superposed do not match.'
1024 do i=0,nsup-(nct-nnt+1)
1025 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1027 nstart_sup=nstart_sup+i
1033 & 'Error - sequences to be superposed do not match.'
1036 if (nsup.eq.0) nsup=nct-nnt
1037 if (nstart_sup.eq.0) nstart_sup=nnt
1038 if (nstart_seq.eq.0) nstart_seq=nnt
1039 if(me.eq.king.or..not.out1file)
1040 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1041 & ' nstart_seq=',nstart_seq
1043 c--- Zscore rms -------
1044 if (nz_start.eq.0) nz_start=nnt
1045 if (nz_end.eq.0 .and. nsup.gt.0) then
1047 else if (nz_end.eq.0) then
1050 if(me.eq.king.or..not.out1file)then
1051 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1052 write (iout,*) 'IZ_SC=',iz_sc
1054 c----------------------
1057 if (.not.pdbref) then
1058 call read_angles(inp,*38)
1060 38 write (iout,'(a)') 'Error reading reference structure.'
1062 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1063 stop 'Error reading reference structure'
1065 39 call chainbuild_extconf
1067 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1074 cref(j,i,kkk)=c(j,i)
1077 call contact(.true.,ncont_ref,icont_ref,co)
1081 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1083 if (constr_dist.gt.0) call read_dist_constr
1084 write (iout,*) "After read_dist_constr nhpb",nhpb
1085 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1087 if(me.eq.king.or..not.out1file)
1088 & write (iout,*) 'Contact order:',co
1090 if(me.eq.king.or..not.out1file)
1091 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1094 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1096 if(me.eq.king.or..not.out1file)
1097 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1098 & icont_ref(1,i),' ',
1099 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1103 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1104 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1105 & modecalc.ne.10) then
1106 C If input structure hasn't been supplied from the PDB file read or generate
1108 if (iranconf.eq.0 .and. .not. extconf) then
1109 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1110 & write (iout,'(a)') 'Initial geometry will be read in.'
1112 read(inp,'(8f10.5)',end=36,err=36)
1113 & ((c(l,k),l=1,3),k=1,nres),
1114 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1115 write (iout,*) "Exit READ_CART"
1116 write (iout,'(8f10.5)')
1117 & ((c(l,k),l=1,3),k=1,nres),
1118 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1119 call int_from_cart1(.true.)
1120 write (iout,*) "Finish INT_TO_CART"
1123 dc(j,i)=c(j,i+1)-c(j,i)
1124 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1128 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1130 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1131 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1137 call read_angles(inp,*36)
1138 call chainbuild_extconf
1141 36 write (iout,'(a)') 'Error reading angle file.'
1143 call mpi_finalize( MPI_COMM_WORLD,IERR )
1145 stop 'Error reading angle file.'
1147 else if (extconf) then
1148 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1149 & write (iout,'(a)') 'Extended chain initial geometry.'
1151 theta(i)=90d0*deg2rad
1154 phi(i)=180d0*deg2rad
1157 alph(i)=110d0*deg2rad
1160 omeg(i)=-120d0*deg2rad
1161 if (itype(i).le.0) omeg(i)=-omeg(i)
1163 call chainbuild_extconf
1165 if(me.eq.king.or..not.out1file)
1166 & write (iout,'(a)') 'Random-generated initial geometry.'
1170 if (me.eq.king .or. fg_rank.eq.0 .and. (
1171 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1175 call gen_rand_conf(itmp,*30)
1177 30 write (iout,*) 'Failed to generate random conformation',
1178 & ', itrial=',itrial
1179 write (*,*) 'Processor:',me,
1180 & ' Failed to generate random conformation',
1190 write (iout,'(a,i3,a)') 'Processor:',me,
1191 & ' error in generating random conformation.'
1192 write (*,'(a,i3,a)') 'Processor:',me,
1193 & ' error in generating random conformation.'
1196 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1201 & ' error in generating random conformation.'
1206 elseif (modecalc.eq.4) then
1207 read (inp,'(a)') intinname
1208 open (intin,file=intinname,status='old',err=333)
1209 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1210 & write (iout,'(a)') 'intinname',intinname
1211 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1213 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1215 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1217 stop 'Error opening angle file.'
1221 C Generate distance constraints, if the PDB structure is to be regularized.
1222 if (nthread.gt.0) then
1223 call read_threadbase
1226 if (me.eq.king .or. .not. out1file)
1228 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1229 write (iout,'(/a,i3,a)')
1230 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1231 write (iout,'(20i4)') (iss(i),i=1,ns)
1233 write(iout,*)"Running with dynamic disulfide-bond formation"
1235 write (iout,'(/a/)') 'Pre-formed links are:'
1241 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1242 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1248 if (ns.gt.0.and.dyn_ss) then
1252 forcon(i-nss)=forcon(i)
1259 dyn_ss_mask(iss(i))=.true.
1262 if (i2ndstr.gt.0) call secstrp2dihc
1263 c call geom_to_var(nvar,x)
1264 c call etotal(energia(0))
1265 c call enerprint(energia(0))
1266 c call briefout(0,etot)
1268 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1269 cd write (iout,'(a)') 'Variable list:'
1270 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1272 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1273 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1274 & 'Processor',myrank,': end reading molecular data.'
1279 c--------------------------------------------------------------------------
1280 logical function seq_comp(itypea,itypeb,length)
1282 integer length,itypea(length),itypeb(length)
1285 if (itypea(i).ne.itypeb(i)) then
1293 c-----------------------------------------------------------------------------
1294 subroutine read_bridge
1295 C Read information about disulfide bridges.
1296 implicit real*8 (a-h,o-z)
1297 include 'DIMENSIONS'
1301 include 'COMMON.IOUNITS'
1302 include 'COMMON.GEO'
1303 include 'COMMON.VAR'
1304 include 'COMMON.INTERACT'
1305 include 'COMMON.LOCAL'
1306 include 'COMMON.NAMES'
1307 include 'COMMON.CHAIN'
1308 include 'COMMON.FFIELD'
1309 include 'COMMON.SBRIDGE'
1310 include 'COMMON.HEADER'
1311 include 'COMMON.CONTROL'
1312 include 'COMMON.DBASE'
1313 include 'COMMON.THREAD'
1314 include 'COMMON.TIME1'
1315 include 'COMMON.SETUP'
1316 C Read bridging residues.
1317 read (inp,*) ns,(iss(i),i=1,ns)
1319 if(me.eq.king.or..not.out1file)
1320 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1321 C Check whether the specified bridging residues are cystines.
1323 if (itype(iss(i)).ne.1) then
1324 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1325 & 'Do you REALLY think that the residue ',
1326 & restyp(itype(iss(i))),i,
1327 & ' can form a disulfide bridge?!!!'
1328 write (*,'(2a,i3,a)')
1329 & 'Do you REALLY think that the residue ',
1330 & restyp(itype(iss(i))),i,
1331 & ' can form a disulfide bridge?!!!'
1333 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1338 C Read preformed bridges.
1340 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1342 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1345 C Check if the residues involved in bridges are in the specified list of
1346 C bridging residues.
1349 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1350 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1351 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1352 & ' contains residues present in other pairs.'
1353 write (*,'(a,i3,a)') 'Disulfide pair',i,
1354 & ' contains residues present in other pairs.'
1356 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1362 if (ihpb(i).eq.iss(j)) goto 10
1364 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1367 if (jhpb(i).eq.iss(j)) goto 20
1369 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1375 ihpb(i)=ihpb(i)+nres
1376 jhpb(i)=jhpb(i)+nres
1382 c----------------------------------------------------------------------------
1383 subroutine read_x(kanal,*)
1384 implicit real*8 (a-h,o-z)
1385 include 'DIMENSIONS'
1386 include 'COMMON.GEO'
1387 include 'COMMON.VAR'
1388 include 'COMMON.CHAIN'
1389 include 'COMMON.IOUNITS'
1390 include 'COMMON.CONTROL'
1391 include 'COMMON.LOCAL'
1392 include 'COMMON.INTERACT'
1393 c Read coordinates from input
1395 read(kanal,'(8f10.5)',end=10,err=10)
1396 & ((c(l,k),l=1,3),k=1,nres),
1397 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1400 c(j,2*nres)=c(j,nres)
1402 call int_from_cart1(.false.)
1405 dc(j,i)=c(j,i+1)-c(j,i)
1406 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1410 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1412 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1413 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1421 c----------------------------------------------------------------------------
1422 subroutine read_threadbase
1423 implicit real*8 (a-h,o-z)
1424 include 'DIMENSIONS'
1425 include 'COMMON.IOUNITS'
1426 include 'COMMON.GEO'
1427 include 'COMMON.VAR'
1428 include 'COMMON.INTERACT'
1429 include 'COMMON.LOCAL'
1430 include 'COMMON.NAMES'
1431 include 'COMMON.CHAIN'
1432 include 'COMMON.FFIELD'
1433 include 'COMMON.SBRIDGE'
1434 include 'COMMON.HEADER'
1435 include 'COMMON.CONTROL'
1436 include 'COMMON.DBASE'
1437 include 'COMMON.THREAD'
1438 include 'COMMON.TIME1'
1439 C Read pattern database for threading.
1440 read (icbase,*) nseq
1442 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1443 & nres_base(2,i),nres_base(3,i)
1444 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1446 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1447 c & nres_base(2,i),nres_base(3,i)
1448 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1452 if (weidis.eq.0.0D0) weidis=0.1D0
1461 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1462 write (iout,'(a,i5)') 'nexcl: ',nexcl
1463 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1466 c------------------------------------------------------------------------------
1467 subroutine setup_var
1468 implicit real*8 (a-h,o-z)
1469 include 'DIMENSIONS'
1470 include 'COMMON.IOUNITS'
1471 include 'COMMON.GEO'
1472 include 'COMMON.VAR'
1473 include 'COMMON.INTERACT'
1474 include 'COMMON.LOCAL'
1475 include 'COMMON.NAMES'
1476 include 'COMMON.CHAIN'
1477 include 'COMMON.FFIELD'
1478 include 'COMMON.SBRIDGE'
1479 include 'COMMON.HEADER'
1480 include 'COMMON.CONTROL'
1481 include 'COMMON.DBASE'
1482 include 'COMMON.THREAD'
1483 include 'COMMON.TIME1'
1484 C Set up variable list.
1490 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1492 ialph(i,1)=nvar+nside
1496 if (indphi.gt.0) then
1498 else if (indback.gt.0) then
1503 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1506 c----------------------------------------------------------------------------
1507 subroutine gen_dist_constr
1508 C Generate CA distance constraints.
1509 implicit real*8 (a-h,o-z)
1510 include 'DIMENSIONS'
1511 include 'COMMON.IOUNITS'
1512 include 'COMMON.GEO'
1513 include 'COMMON.VAR'
1514 include 'COMMON.INTERACT'
1515 include 'COMMON.LOCAL'
1516 include 'COMMON.NAMES'
1517 include 'COMMON.CHAIN'
1518 include 'COMMON.FFIELD'
1519 include 'COMMON.SBRIDGE'
1520 include 'COMMON.HEADER'
1521 include 'COMMON.CONTROL'
1522 include 'COMMON.DBASE'
1523 include 'COMMON.THREAD'
1524 include 'COMMON.TIME1'
1525 dimension itype_pdb(maxres)
1526 common /pizda/ itype_pdb
1528 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1529 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1530 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1532 if (constr_dist.eq.11) then
1533 do i=nstart_sup,nstart_sup+nsup-1
1534 do j=i+2,nstart_sup+nsup-1
1536 if (distance.le.15.0) then
1538 ihpb(nhpb)=i+nstart_seq-nstart_sup
1539 jhpb(nhpb)=j+nstart_seq-nstart_sup
1540 forcon(nhpb)=sqrt(0.04*distance)
1541 fordepth(nhpb)=sqrt(40.0/distance)
1542 dhpb(nhpb)=distance-0.1d0
1543 dhpb1(nhpb)=distance+0.1d0
1546 if (.not.out1file .or. me.eq.king)
1547 & write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1548 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1550 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1551 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1557 do i=nstart_sup,nstart_sup+nsup-1
1558 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1559 cd & ' seq_pdb', restyp(itype_pdb(i))
1560 do j=i+2,nstart_sup+nsup-1
1562 ihpb(nhpb)=i+nstart_seq-nstart_sup
1563 jhpb(nhpb)=j+nstart_seq-nstart_sup
1565 dhpb(nhpb)=dist(i,j)
1569 cd write (iout,'(a)') 'Distance constraints:'
1574 cd if (ii.gt.nres) then
1579 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1580 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1581 cd & dhpb(i),forcon(i)
1585 c----------------------------------------------------------------------------
1587 implicit real*8 (a-h,o-z)
1588 include 'DIMENSIONS'
1589 include 'COMMON.MAP'
1590 include 'COMMON.IOUNITS'
1591 character*3 angid(4) /'THE','PHI','ALP','OME'/
1592 character*80 mapcard,ucase
1594 read (inp,'(a)') mapcard
1595 mapcard=ucase(mapcard)
1596 if (index(mapcard,'PHI').gt.0) then
1598 else if (index(mapcard,'THE').gt.0) then
1600 else if (index(mapcard,'ALP').gt.0) then
1602 else if (index(mapcard,'OME').gt.0) then
1605 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1606 stop 'Error - illegal variable spec in MAP card.'
1608 call readi (mapcard,'RES1',res1(imap),0)
1609 call readi (mapcard,'RES2',res2(imap),0)
1610 if (res1(imap).eq.0) then
1611 res1(imap)=res2(imap)
1612 else if (res2(imap).eq.0) then
1613 res2(imap)=res1(imap)
1615 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1617 & 'Error - illegal definition of variable group in MAP.'
1618 stop 'Error - illegal definition of variable group in MAP.'
1620 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1621 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1622 call readi(mapcard,'NSTEP',nstep(imap),0)
1623 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1625 & 'Illegal boundary and/or step size specification in MAP.'
1626 stop 'Illegal boundary and/or step size specification in MAP.'
1631 c----------------------------------------------------------------------------
1633 implicit real*8 (a-h,o-z)
1634 include 'DIMENSIONS'
1635 include 'COMMON.IOUNITS'
1636 include 'COMMON.GEO'
1637 include 'COMMON.CSA'
1638 include 'COMMON.BANK'
1639 include 'COMMON.CONTROL'
1641 character*620 mcmcard
1642 call card_concat(mcmcard)
1644 call readi(mcmcard,'NCONF',nconf,50)
1645 call readi(mcmcard,'NADD',nadd,0)
1646 call readi(mcmcard,'JSTART',jstart,1)
1647 call readi(mcmcard,'JEND',jend,1)
1648 call readi(mcmcard,'NSTMAX',nstmax,500000)
1649 call readi(mcmcard,'N0',n0,1)
1650 call readi(mcmcard,'N1',n1,6)
1651 call readi(mcmcard,'N2',n2,4)
1652 call readi(mcmcard,'N3',n3,0)
1653 call readi(mcmcard,'N4',n4,0)
1654 call readi(mcmcard,'N5',n5,0)
1655 call readi(mcmcard,'N6',n6,10)
1656 call readi(mcmcard,'N7',n7,0)
1657 call readi(mcmcard,'N8',n8,0)
1658 call readi(mcmcard,'N9',n9,0)
1659 call readi(mcmcard,'N14',n14,0)
1660 call readi(mcmcard,'N15',n15,0)
1661 call readi(mcmcard,'N16',n16,0)
1662 call readi(mcmcard,'N17',n17,0)
1663 call readi(mcmcard,'N18',n18,0)
1665 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1667 call readi(mcmcard,'NDIFF',ndiff,2)
1668 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1669 call readi(mcmcard,'IS1',is1,1)
1670 call readi(mcmcard,'IS2',is2,8)
1671 call readi(mcmcard,'NRAN0',nran0,4)
1672 call readi(mcmcard,'NRAN1',nran1,2)
1673 call readi(mcmcard,'IRR',irr,1)
1674 call readi(mcmcard,'NSEED',nseed,20)
1675 call readi(mcmcard,'NTOTAL',ntotal,10000)
1676 call reada(mcmcard,'CUT1',cut1,2.0d0)
1677 call reada(mcmcard,'CUT2',cut2,5.0d0)
1678 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1679 call readi(mcmcard,'ICMAX',icmax,3)
1680 call readi(mcmcard,'IRESTART',irestart,0)
1681 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1684 call reada(mcmcard,'DELE',dele,20.0d0)
1685 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1686 call readi(mcmcard,'IREF',iref,0)
1687 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1688 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1689 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1690 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1691 write (iout,*) "NCONF_IN",nconf_in
1694 c----------------------------------------------------------------------------
1695 cfmc subroutine mcmfread
1696 cfmc implicit real*8 (a-h,o-z)
1697 cfmc include 'DIMENSIONS'
1698 cfmc include 'COMMON.MCMF'
1699 cfmc include 'COMMON.IOUNITS'
1700 cfmc include 'COMMON.GEO'
1701 cfmc character*80 ucase
1702 cfmc character*620 mcmcard
1703 cfmc call card_concat(mcmcard)
1705 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1706 cfmc write(iout,*)'MAXRANT=',maxrant
1707 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1708 cfmc write(iout,*)'MAXFAM=',maxfam
1709 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1710 cfmc write(iout,*)'NNET1=',nnet1
1711 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1712 cfmc write(iout,*)'NNET2=',nnet2
1713 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1714 cfmc write(iout,*)'NNET3=',nnet3
1715 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1716 cfmc write(iout,*)'ILASTT=',ilastt
1717 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1718 cfmc write(iout,*)'MAXSTR=',maxstr
1719 cfmc maxstr_f=maxstr/maxfam
1720 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1721 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1722 cfmc write(iout,*)'NMCMF=',nmcmf
1723 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1724 cfmc write(iout,*)'IFOCUS=',ifocus
1725 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1726 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1727 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1728 cfmc write(iout,*)'INTPRT=',intprt
1729 cfmc call readi(mcmcard,'IPRT',iprt,100)
1730 cfmc write(iout,*)'IPRT=',iprt
1731 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1732 cfmc write(iout,*)'IMAXTR=',imaxtr
1733 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1734 cfmc write(iout,*)'MAXEVEN=',maxeven
1735 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1736 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1737 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1738 cfmc write(iout,*)'INIMIN=',inimin
1739 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1740 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1741 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1742 cfmc write(iout,*)'NTHREAD=',nthread
1743 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1744 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1745 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1746 cfmc write(iout,*)'MAXPERT=',maxpert
1747 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1748 cfmc write(iout,*)'IRMSD=',irmsd
1749 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1750 cfmc write(iout,*)'DENEMIN=',denemin
1751 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1752 cfmc write(iout,*)'RCUT1S=',rcut1s
1753 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1754 cfmc write(iout,*)'RCUT1E=',rcut1e
1755 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1756 cfmc write(iout,*)'RCUT2S=',rcut2s
1757 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1758 cfmc write(iout,*)'RCUT2E=',rcut2e
1759 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1760 cfmc write(iout,*)'DPERT1=',d_pert1
1761 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1762 cfmc write(iout,*)'DPERT1A=',d_pert1a
1763 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1764 cfmc write(iout,*)'DPERT2=',d_pert2
1765 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1766 cfmc write(iout,*)'DPERT2A=',d_pert2a
1767 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1768 cfmc write(iout,*)'DPERT2B=',d_pert2b
1769 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1770 cfmc write(iout,*)'DPERT2C=',d_pert2c
1771 cfmc d_pert1=deg2rad*d_pert1
1772 cfmc d_pert1a=deg2rad*d_pert1a
1773 cfmc d_pert2=deg2rad*d_pert2
1774 cfmc d_pert2a=deg2rad*d_pert2a
1775 cfmc d_pert2b=deg2rad*d_pert2b
1776 cfmc d_pert2c=deg2rad*d_pert2c
1777 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1778 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1779 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1780 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1781 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1782 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1783 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1784 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1785 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1786 cfmc write(iout,*)'RCUTINI=',rcutini
1787 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1788 cfmc write(iout,*)'GRAT=',grat
1789 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1790 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1794 c----------------------------------------------------------------------------
1796 implicit real*8 (a-h,o-z)
1797 include 'DIMENSIONS'
1798 include 'COMMON.MCM'
1799 include 'COMMON.MCE'
1800 include 'COMMON.IOUNITS'
1802 character*320 mcmcard
1803 call card_concat(mcmcard)
1804 call readi(mcmcard,'MAXACC',maxacc,100)
1805 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1806 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1807 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1808 call readi(mcmcard,'MAXREPM',maxrepm,200)
1809 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1810 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1811 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1812 call reada(mcmcard,'E_UP',e_up,5.0D0)
1813 call reada(mcmcard,'DELTE',delte,0.1D0)
1814 call readi(mcmcard,'NSWEEP',nsweep,5)
1815 call readi(mcmcard,'NSTEPH',nsteph,0)
1816 call readi(mcmcard,'NSTEPC',nstepc,0)
1817 call reada(mcmcard,'TMIN',tmin,298.0D0)
1818 call reada(mcmcard,'TMAX',tmax,298.0D0)
1819 call readi(mcmcard,'NWINDOW',nwindow,0)
1820 call readi(mcmcard,'PRINT_MC',print_mc,0)
1821 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1822 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1823 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1824 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1825 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1826 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1827 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1828 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1829 if (nwindow.gt.0) then
1830 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1832 winlen(i)=winend(i)-winstart(i)+1
1835 if (tmax.lt.tmin) tmax=tmin
1836 if (tmax.eq.tmin) then
1840 if (nstepc.gt.0 .and. nsteph.gt.0) then
1841 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1842 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1844 C Probabilities of different move types
1845 sumpro_type(0)=0.0D0
1846 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1847 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1848 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1849 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1850 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1851 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1852 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1854 print *,'i',i,' sumprotype',sumpro_type(i)
1855 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1856 print *,'i',i,' sumprotype',sumpro_type(i)
1860 c----------------------------------------------------------------------------
1861 subroutine read_minim
1862 implicit real*8 (a-h,o-z)
1863 include 'DIMENSIONS'
1864 include 'COMMON.MINIM'
1865 include 'COMMON.IOUNITS'
1867 character*320 minimcard
1868 call card_concat(minimcard)
1869 call readi(minimcard,'MAXMIN',maxmin,2000)
1870 call readi(minimcard,'MAXFUN',maxfun,5000)
1871 call readi(minimcard,'MINMIN',minmin,maxmin)
1872 call readi(minimcard,'MINFUN',minfun,maxmin)
1873 call reada(minimcard,'TOLF',tolf,1.0D-2)
1874 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1875 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1876 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1877 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1878 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1879 & 'Options in energy minimization:'
1880 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1881 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1882 & 'MinMin:',MinMin,' MinFun:',MinFun,
1883 & ' TolF:',TolF,' RTolF:',RTolF
1886 c----------------------------------------------------------------------------
1887 subroutine read_angles(kanal,*)
1888 implicit real*8 (a-h,o-z)
1889 include 'DIMENSIONS'
1890 include 'COMMON.GEO'
1891 include 'COMMON.VAR'
1892 include 'COMMON.CHAIN'
1893 include 'COMMON.IOUNITS'
1894 include 'COMMON.CONTROL'
1895 c Read angles from input
1897 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1898 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1899 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1900 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1903 c 9/7/01 avoid 180 deg valence angle
1904 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1906 theta(i)=deg2rad*theta(i)
1907 phi(i)=deg2rad*phi(i)
1908 alph(i)=deg2rad*alph(i)
1909 omeg(i)=deg2rad*omeg(i)
1914 c----------------------------------------------------------------------------
1915 subroutine reada(rekord,lancuch,wartosc,default)
1917 character*(*) rekord,lancuch
1918 double precision wartosc,default
1921 iread=index(rekord,lancuch)
1922 if (iread.eq.0) then
1926 iread=iread+ilen(lancuch)+1
1927 read (rekord(iread:),*,err=10,end=10) wartosc
1932 c----------------------------------------------------------------------------
1933 subroutine readi(rekord,lancuch,wartosc,default)
1935 character*(*) rekord,lancuch
1936 integer wartosc,default
1939 iread=index(rekord,lancuch)
1940 if (iread.eq.0) then
1944 iread=iread+ilen(lancuch)+1
1945 read (rekord(iread:),*,err=10,end=10) wartosc
1950 c----------------------------------------------------------------------------
1951 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1954 integer tablica(dim),default
1955 character*(*) rekord,lancuch
1962 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1963 if (iread.eq.0) return
1964 iread=iread+ilen(lancuch)+1
1965 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1968 c----------------------------------------------------------------------------
1969 subroutine multreada(rekord,lancuch,tablica,dim,default)
1972 double precision tablica(dim),default
1973 character*(*) rekord,lancuch
1980 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1981 if (iread.eq.0) return
1982 iread=iread+ilen(lancuch)+1
1983 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1986 c----------------------------------------------------------------------------
1987 subroutine openunits
1988 implicit real*8 (a-h,o-z)
1989 include 'DIMENSIONS'
1992 character*16 form,nodename
1995 include 'COMMON.SETUP'
1996 include 'COMMON.IOUNITS'
1998 include 'COMMON.CONTROL'
1999 integer lenpre,lenpot,ilen,lentmp
2001 character*3 out1file_text,ucase
2004 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2005 call getenv_loc("PREFIX",prefix)
2007 call getenv_loc("POT",pot)
2008 call getenv_loc("DIRTMP",tmpdir)
2009 call getenv_loc("CURDIR",curdir)
2010 call getenv_loc("OUT1FILE",out1file_text)
2011 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2012 out1file_text=ucase(out1file_text)
2013 if (out1file_text(1:1).eq."Y") then
2016 out1file=fg_rank.gt.0
2021 if (lentmp.gt.0) then
2022 write (*,'(80(1h!))')
2023 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2024 write (*,'(80(1h!))')
2025 write (*,*)"All output files will be on node /tmp directory."
2027 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2028 if (me.eq.king) then
2029 write (*,*) "The master node is ",nodename
2030 else if (fg_rank.eq.0) then
2031 write (*,*) "I am the CG slave node ",nodename
2033 write (*,*) "I am the FG slave node ",nodename
2036 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2037 lenpre = lentmp+lenpre+1
2039 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2040 C Get the names and open the input files
2041 #if defined(WINIFL) || defined(WINPGI)
2042 open(1,file=pref_orig(:ilen(pref_orig))//
2043 & '.inp',status='old',readonly,shared)
2044 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2045 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2046 C Get parameter filenames and open the parameter files.
2047 call getenv_loc('BONDPAR',bondname)
2048 open (ibond,file=bondname,status='old',readonly,shared)
2049 call getenv_loc('THETPAR',thetname)
2050 open (ithep,file=thetname,status='old',readonly,shared)
2051 call getenv_loc('ROTPAR',rotname)
2052 open (irotam,file=rotname,status='old',readonly,shared)
2053 call getenv_loc('TORPAR',torname)
2054 open (itorp,file=torname,status='old',readonly,shared)
2055 call getenv_loc('TORDPAR',tordname)
2056 open (itordp,file=tordname,status='old',readonly,shared)
2057 call getenv_loc('TORKCC',torkccname)
2058 open (itorkcc,file=torkccname,status='old',readonly,shared)
2059 call getenv_loc('THETKCC',thetkccname)
2060 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2061 call getenv_loc('FOURIER',fouriername)
2062 open (ifourier,file=fouriername,status='old',readonly,shared)
2063 call getenv_loc('ELEPAR',elename)
2064 open (ielep,file=elename,status='old',readonly,shared)
2065 call getenv_loc('SIDEPAR',sidename)
2066 open (isidep,file=sidename,status='old',readonly,shared)
2067 call getenv_loc('LIPTRANPAR',liptranname)
2068 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2069 #elif (defined CRAY) || (defined AIX)
2070 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2072 c print *,"Processor",myrank," opened file 1"
2073 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2074 c print *,"Processor",myrank," opened file 9"
2075 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2076 C Get parameter filenames and open the parameter files.
2077 call getenv_loc('BONDPAR',bondname)
2078 open (ibond,file=bondname,status='old',action='read')
2079 c print *,"Processor",myrank," opened file IBOND"
2080 call getenv_loc('THETPAR',thetname)
2081 open (ithep,file=thetname,status='old',action='read')
2082 c print *,"Processor",myrank," opened file ITHEP"
2083 call getenv_loc('ROTPAR',rotname)
2084 open (irotam,file=rotname,status='old',action='read')
2085 c print *,"Processor",myrank," opened file IROTAM"
2086 call getenv_loc('TORPAR',torname)
2087 open (itorp,file=torname,status='old',action='read')
2088 c print *,"Processor",myrank," opened file ITORP"
2089 call getenv_loc('TORDPAR',tordname)
2090 open (itordp,file=tordname,status='old',action='read')
2091 call getenv_loc('TORKCC',torkccname)
2092 open (itorkcc,file=torkccname,status='old',action='read')
2093 call getenv_loc('THETKCC',thetkccname)
2094 open (ithetkcc,file=thetkccname,status='old',action='read')
2095 c print *,"Processor",myrank," opened file ITORDP"
2096 call getenv_loc('SCCORPAR',sccorname)
2097 open (isccor,file=sccorname,status='old',action='read')
2098 c print *,"Processor",myrank," opened file ISCCOR"
2099 call getenv_loc('FOURIER',fouriername)
2100 open (ifourier,file=fouriername,status='old',action='read')
2101 c print *,"Processor",myrank," opened file IFOURIER"
2102 call getenv_loc('ELEPAR',elename)
2103 open (ielep,file=elename,status='old',action='read')
2104 c print *,"Processor",myrank," opened file IELEP"
2105 call getenv_loc('SIDEPAR',sidename)
2106 open (isidep,file=sidename,status='old',action='read')
2107 call getenv_loc('LIPTRANPAR',liptranname)
2108 open (iliptranpar,file=liptranname,status='old',action='read')
2109 c print *,"Processor",myrank," opened file ISIDEP"
2110 c print *,"Processor",myrank," opened parameter files"
2112 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2113 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2114 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2115 C Get parameter filenames and open the parameter files.
2116 call getenv_loc('BONDPAR',bondname)
2117 open (ibond,file=bondname,status='old')
2118 call getenv_loc('THETPAR',thetname)
2119 open (ithep,file=thetname,status='old')
2120 call getenv_loc('ROTPAR',rotname)
2121 open (irotam,file=rotname,status='old')
2122 call getenv_loc('TORPAR',torname)
2123 open (itorp,file=torname,status='old')
2124 call getenv_loc('TORDPAR',tordname)
2125 open (itordp,file=tordname,status='old')
2126 call getenv_loc('TORKCC',torkccname)
2127 open (itorkcc,file=torkccname,status='old')
2128 call getenv_loc('THETKCC',thetkccname)
2129 open (ithetkcc,file=thetkccname,status='old')
2130 call getenv_loc('SCCORPAR',sccorname)
2131 open (isccor,file=sccorname,status='old')
2132 call getenv_loc('FOURIER',fouriername)
2133 open (ifourier,file=fouriername,status='old')
2134 call getenv_loc('ELEPAR',elename)
2135 open (ielep,file=elename,status='old')
2136 call getenv_loc('SIDEPAR',sidename)
2137 open (isidep,file=sidename,status='old')
2138 call getenv_loc('LIPTRANPAR',liptranname)
2139 open (iliptranpar,file=liptranname,status='old')
2141 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2143 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2144 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2145 C Get parameter filenames and open the parameter files.
2146 call getenv_loc('BONDPAR',bondname)
2147 open (ibond,file=bondname,status='old',readonly)
2148 call getenv_loc('THETPAR',thetname)
2149 open (ithep,file=thetname,status='old',readonly)
2150 call getenv_loc('ROTPAR',rotname)
2151 open (irotam,file=rotname,status='old',readonly)
2152 call getenv_loc('TORPAR',torname)
2153 open (itorp,file=torname,status='old',readonly)
2154 call getenv_loc('TORDPAR',tordname)
2155 open (itordp,file=tordname,status='old',readonly)
2156 call getenv_loc('TORKCC',torkccname)
2157 open (itorkcc,file=torkccname,status='old',readonly)
2158 call getenv_loc('THETKCC',thetkccname)
2159 open (ithetkcc,file=thetkccname,status='old',readonly)
2160 call getenv_loc('SCCORPAR',sccorname)
2161 open (isccor,file=sccorname,status='old',readonly)
2163 call getenv_loc('THETPARPDB',thetname_pdb)
2164 print *,"thetname_pdb ",thetname_pdb
2165 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2166 print *,ithep_pdb," opened"
2168 call getenv_loc('FOURIER',fouriername)
2169 open (ifourier,file=fouriername,status='old',readonly)
2170 call getenv_loc('ELEPAR',elename)
2171 open (ielep,file=elename,status='old',readonly)
2172 call getenv_loc('SIDEPAR',sidename)
2173 open (isidep,file=sidename,status='old',readonly)
2174 call getenv_loc('LIPTRANPAR',liptranname)
2175 open (iliptranpar,file=liptranname,status='old',action='read')
2177 call getenv_loc('ROTPARPDB',rotname_pdb)
2178 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2181 call getenv_loc('TUBEPAR',tubename)
2182 open (itube,file=tubename,status='old',readonly)
2185 C 8/9/01 In the newest version SCp interaction constants are read from a file
2186 C Use -DOLDSCP to use hard-coded constants instead.
2188 call getenv_loc('SCPPAR',scpname)
2189 #if defined(WINIFL) || defined(WINPGI)
2190 open (iscpp,file=scpname,status='old',readonly,shared)
2191 #elif (defined CRAY) || (defined AIX)
2192 open (iscpp,file=scpname,status='old',action='read')
2194 open (iscpp,file=scpname,status='old')
2196 open (iscpp,file=scpname,status='old',readonly)
2199 call getenv_loc('PATTERN',patname)
2200 #if defined(WINIFL) || defined(WINPGI)
2201 open (icbase,file=patname,status='old',readonly,shared)
2202 #elif (defined CRAY) || (defined AIX)
2203 open (icbase,file=patname,status='old',action='read')
2205 open (icbase,file=patname,status='old')
2207 open (icbase,file=patname,status='old',readonly)
2210 C Open output file only for CG processes
2211 c print *,"Processor",myrank," fg_rank",fg_rank
2212 if (fg_rank.eq.0) then
2214 if (nodes.eq.1) then
2217 npos = dlog10(dfloat(nodes-1))+1
2219 if (npos.lt.3) npos=3
2220 write (liczba,'(i1)') npos
2221 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2223 write (liczba,form) me
2224 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2225 & liczba(:ilen(liczba))
2226 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2228 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2230 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2231 & liczba(:ilen(liczba))//'.mol2'
2232 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2233 & liczba(:ilen(liczba))//'.stat'
2235 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2236 & //liczba(:ilen(liczba))//'.stat')
2237 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2240 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2241 & liczba(:ilen(liczba))//'.const'
2246 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2247 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2248 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2249 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2250 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2252 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2254 rest2name=prefix(:ilen(prefix))//'.rst'
2256 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2259 #if defined(AIX) || defined(PGI)
2260 if (me.eq.king .or. .not. out1file)
2261 & open(iout,file=outname,status='unknown')
2263 if (fg_rank.gt.0) then
2264 write (liczba,'(i3.3)') myrank/nfgtasks
2265 write (ll,'(bz,i3.3)') fg_rank
2266 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2271 open(igeom,file=intname,status='unknown',position='append')
2272 open(ipdb,file=pdbname,status='unknown')
2273 open(imol2,file=mol2name,status='unknown')
2274 open(istat,file=statname,status='unknown',position='append')
2276 c1out open(iout,file=outname,status='unknown')
2279 if (me.eq.king .or. .not.out1file)
2280 & open(iout,file=outname,status='unknown')
2282 if (fg_rank.gt.0) then
2283 write (liczba,'(i3.3)') myrank/nfgtasks
2284 write (ll,'(bz,i3.3)') fg_rank
2285 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2290 open(igeom,file=intname,status='unknown',access='append')
2291 open(ipdb,file=pdbname,status='unknown')
2292 open(imol2,file=mol2name,status='unknown')
2293 open(istat,file=statname,status='unknown',access='append')
2295 c1out open(iout,file=outname,status='unknown')
2298 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2299 csa_seed=prefix(:lenpre)//'.CSA.seed'
2300 csa_history=prefix(:lenpre)//'.CSA.history'
2301 csa_bank=prefix(:lenpre)//'.CSA.bank'
2302 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2303 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2304 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2305 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2306 csa_int=prefix(:lenpre)//'.int'
2307 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2308 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2309 csa_in=prefix(:lenpre)//'.CSA.in'
2310 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2313 write (iout,'(80(1h-))')
2314 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2315 write (iout,'(80(1h-))')
2316 write (iout,*) "Input file : ",
2317 & pref_orig(:ilen(pref_orig))//'.inp'
2318 write (iout,*) "Output file : ",
2319 & outname(:ilen(outname))
2321 write (iout,*) "Sidechain potential file : ",
2322 & sidename(:ilen(sidename))
2324 write (iout,*) "SCp potential file : ",
2325 & scpname(:ilen(scpname))
2327 write (iout,*) "Electrostatic potential file : ",
2328 & elename(:ilen(elename))
2329 write (iout,*) "Cumulant coefficient file : ",
2330 & fouriername(:ilen(fouriername))
2331 write (iout,*) "Torsional parameter file : ",
2332 & torname(:ilen(torname))
2333 write (iout,*) "Double torsional parameter file : ",
2334 & tordname(:ilen(tordname))
2335 write (iout,*) "SCCOR parameter file : ",
2336 & sccorname(:ilen(sccorname))
2337 write (iout,*) "Bond & inertia constant file : ",
2338 & bondname(:ilen(bondname))
2339 write (iout,*) "Bending parameter file : ",
2340 & thetname(:ilen(thetname))
2341 write (iout,*) "Rotamer parameter file : ",
2342 & rotname(:ilen(rotname))
2343 write (iout,*) "Thetpdb parameter file : ",
2344 & thetname_pdb(:ilen(thetname_pdb))
2345 write (iout,*) "Threading database : ",
2346 & patname(:ilen(patname))
2348 &write (iout,*)" DIRTMP : ",
2350 write (iout,'(80(1h-))')
2354 c----------------------------------------------------------------------------
2355 subroutine card_concat(card)
2356 implicit real*8 (a-h,o-z)
2357 include 'DIMENSIONS'
2358 include 'COMMON.IOUNITS'
2360 character*80 karta,ucase
2362 read (inp,'(a)') karta
2365 do while (karta(80:80).eq.'&')
2366 card=card(:ilen(card)+1)//karta(:79)
2367 read (inp,'(a)') karta
2370 card=card(:ilen(card)+1)//karta
2373 c----------------------------------------------------------------------------------
2375 implicit real*8 (a-h,o-z)
2376 include 'DIMENSIONS'
2377 include 'COMMON.CHAIN'
2378 include 'COMMON.IOUNITS'
2380 open(irest2,file=rest2name,status='unknown')
2381 read(irest2,*) totT,EK,potE,totE,t_bath
2384 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2387 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2390 read (irest2,*) iset
2395 c---------------------------------------------------------------------------------
2396 subroutine read_fragments
2397 implicit real*8 (a-h,o-z)
2398 include 'DIMENSIONS'
2402 include 'COMMON.SETUP'
2403 include 'COMMON.CHAIN'
2404 include 'COMMON.IOUNITS'
2406 include 'COMMON.CONTROL'
2407 read(inp,*) nset,nfrag,npair,nfrag_back
2408 if(me.eq.king.or..not.out1file)
2409 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2410 & " nfrag_back",nfrag_back
2412 read(inp,*) mset(iset)
2414 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2416 if(me.eq.king.or..not.out1file)
2417 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2418 & ifrag(2,i,iset), qinfrag(i,iset)
2421 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2423 if(me.eq.king.or..not.out1file)
2424 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2425 & ipair(2,i,iset), qinpair(i,iset)
2428 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2429 & wfrag_back(3,i,iset),
2430 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2431 if(me.eq.king.or..not.out1file)
2432 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2433 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2438 C---------------------------------------------------------------------------
2439 subroutine read_afminp
2440 implicit real*8 (a-h,o-z)
2441 include 'DIMENSIONS'
2445 include 'COMMON.SETUP'
2446 include 'COMMON.CONTROL'
2447 include 'COMMON.CHAIN'
2448 include 'COMMON.IOUNITS'
2449 include 'COMMON.SBRIDGE'
2450 character*320 afmcard
2452 call card_concat(afmcard)
2453 call readi(afmcard,"BEG",afmbeg,0)
2454 call readi(afmcard,"END",afmend,0)
2455 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2456 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2457 print *,'FORCE=' ,forceAFMconst
2458 CCCC NOW PROPERTIES FOR AFM
2461 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2463 distafminit=dsqrt(distafminit)
2464 print *,'initdist',distafminit
2468 c-------------------------------------------------------------------------------
2469 subroutine read_dist_constr
2470 implicit real*8 (a-h,o-z)
2471 include 'DIMENSIONS'
2475 include 'COMMON.SETUP'
2476 include 'COMMON.CONTROL'
2477 include 'COMMON.CHAIN'
2478 include 'COMMON.IOUNITS'
2479 include 'COMMON.SBRIDGE'
2480 integer ifrag_(2,100),ipair_(2,100)
2481 double precision wfrag_(100),wpair_(100)
2482 character*500 controlcard
2484 write (iout,*) "Calling read_dist_constr"
2485 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2487 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
2488 call gen_dist_constr
2491 call card_concat(controlcard)
2492 call readi(controlcard,"NFRAG",nfrag_,0)
2493 call readi(controlcard,"NPAIR",npair_,0)
2494 call readi(controlcard,"NDIST",ndist_,0)
2495 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2496 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2497 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2498 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2499 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2500 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2501 c write (iout,*) "IFRAG"
2503 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2505 c write (iout,*) "IPAIR"
2507 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2511 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2512 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2513 & ifrag_(2,i)=nstart_sup+nsup-1
2514 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2516 if (wfrag_(i).gt.0.0d0) then
2517 do j=ifrag_(1,i),ifrag_(2,i)-1
2518 do k=j+1,ifrag_(2,i)
2519 c write (iout,*) "j",j," k",k
2521 if (constr_dist.eq.1) then
2526 forcon(nhpb)=wfrag_(i)
2527 else if (constr_dist.eq.2) then
2528 if (ddjk.le.dist_cut) then
2533 forcon(nhpb)=wfrag_(i)
2540 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2543 if (.not.out1file .or. me.eq.king)
2544 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2545 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2547 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2548 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2555 if (wpair_(i).gt.0.0d0) then
2563 do j=ifrag_(1,ii),ifrag_(2,ii)
2564 do k=ifrag_(1,jj),ifrag_(2,jj)
2568 forcon(nhpb)=wpair_(i)
2569 dhpb(nhpb)=dist(j,k)
2571 if (.not.out1file .or. me.eq.king)
2572 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2573 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2575 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2576 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2584 if (constr_dist.eq.11) then
2585 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2586 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2587 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2590 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2591 & ibecarb(i),forcon(nhpb+1)
2593 if (forcon(nhpb+1).gt.0.0d0) then
2595 if (ibecarb(i).gt.0) then
2596 ihpb(i)=ihpb(i)+nres
2597 jhpb(i)=jhpb(i)+nres
2599 if (dhpb(nhpb).eq.0.0d0)
2600 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2602 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2603 C if (forcon(nhpb+1).gt.0.0d0) then
2605 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2607 if (.not.out1file .or. me.eq.king)
2608 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2609 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2611 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2612 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2620 c-------------------------------------------------------------------------------
2622 subroutine flush(iu)
2627 subroutine flush(iu)
2632 c------------------------------------------------------------------------------
2633 subroutine copy_to_tmp(source)
2634 include "DIMENSIONS"
2635 include "COMMON.IOUNITS"
2636 character*(*) source
2637 character* 256 tmpfile
2641 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2642 inquire(file=tmpfile,exist=ex)
2644 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2645 & " to temporary directory..."
2646 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2647 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2651 c------------------------------------------------------------------------------
2652 subroutine move_from_tmp(source)
2653 include "DIMENSIONS"
2654 include "COMMON.IOUNITS"
2655 character*(*) source
2658 write (*,*) "Moving ",source(:ilen(source)),
2659 & " from temporary directory to working directory"
2660 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2661 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2664 c------------------------------------------------------------------------------
2665 subroutine random_init(seed)
2667 C Initialize random number generator
2669 implicit real*8 (a-h,o-z)
2670 include 'DIMENSIONS'
2673 logical OKRandom, prng_restart
2675 integer iseed_array(4)
2677 include 'COMMON.IOUNITS'
2678 include 'COMMON.TIME1'
2679 include 'COMMON.THREAD'
2680 include 'COMMON.SBRIDGE'
2681 include 'COMMON.CONTROL'
2682 include 'COMMON.MCM'
2683 include 'COMMON.MAP'
2684 include 'COMMON.HEADER'
2685 include 'COMMON.CSA'
2686 include 'COMMON.CHAIN'
2687 include 'COMMON.MUCA'
2689 include 'COMMON.FFIELD'
2690 include 'COMMON.SETUP'
2691 iseed=-dint(dabs(seed))
2692 if (iseed.eq.0) then
2693 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2694 & 'Random seed undefined. The program will stop.'
2695 write (*,'(/80(1h*)/20x,a/80(1h*))')
2696 & 'Random seed undefined. The program will stop.'
2698 call mpi_finalize(mpi_comm_world,ierr)
2700 stop 'Bad random seed.'
2703 if (fg_rank.eq.0) then
2707 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2708 OKRandom = prng_restart(me,iseed)
2711 tmp=65536.0d0**(4-i)
2712 iseed_array(i) = dint(seed/tmp)
2713 seed=seed-iseed_array(i)*tmp
2716 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2717 & (iseed_array(i),i=1,4)
2718 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2719 & (iseed_array(i),i=1,4)
2720 OKRandom = prng_restart(me,iseed_array)
2723 c r1 = prng_next(me)
2724 r1=ran_number(0.0D0,1.0D0)
2726 & write (iout,*) 'ran_num',r1
2727 if (r1.lt.0.0d0) OKRandom=.false.
2729 if (.not.OKRandom) then
2730 write (iout,*) 'PRNG IS NOT WORKING!!!'
2731 print *,'PRNG IS NOT WORKING!!!'
2734 call mpi_abort(mpi_comm_world,error_msg,ierr)
2737 write (iout,*) 'too many processors for parallel prng'
2738 write (*,*) 'too many processors for parallel prng'
2746 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)